function [ validation_report ] = validate_sicd( sicd_filename, do_interactive )
%VALIDATE_SICD Generates a validation report for a given SICD
%
% This is the master function that orchestrates all the steps of automated
% SICD validation, as well as some user-interactive ones.
%
% It is recommended that when examining a new source of SICD data, all
% modes from that source be reviewed as well as left- and right-looking
% datasets, if available.
%
% This code is currently just a freeform list of all the SICD metadata bugs
% we have discovered browsing through many, many attempts at SICD from
% multiple SAR sources.  This list is certainly not complete, and possibly
% not even all that well organized, but hopefully it serves as a base to
% which other test can be added as they are discovered.
%
% Possible improvements that could be made:
% 1) If the user doesn't want to do full image analysis, we could
% probably make a measurement of the weighting functions by sampling lines
% and rows rather than processing every single pixel in the entire image.
% 2) Spatial frequency to data comparison is current done through manual
% interactive review.  Probably a way to at least partially automate this.
%
% Written by: Wade Schwartzkopf, NGA/IDT
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

%% Setup
validation_report = cell(0,3);
if ~exist('do_interactive','var')
    do_interactive = true;
end
% Extract SICD XML data from file
NITF_meta = read_sicd_nitf_offsets( sicd_filename );
fid = fopen(sicd_filename);
fseek(fid,NITF_meta.minimal.desOffsets(1),'bof');
SICD_xml_string = fread(fid,NITF_meta.minimal.desLengths(1),'uint8=>char')';
fclose(fid);
% Convert metadata into MATLAB structure
SICD_meta = sicdxml2struct( xmlread( java.io.StringBufferInputStream( SICD_xml_string ) ) );

%% 1. Most fundamental check.  First validate XML schema
% SICD XML should validate with one (and only one) schema.  If it does, we
% will report which schema it validated against.  If it does not, we will
% report the first validation error against each schema.
% Now perform validation
factory = javax.xml.validation.SchemaFactory.newInstance('http://www.w3.org/2001/XMLSchema');
xsd_path = fileparts(which(mfilename('fullpath')));
filelist = dir(fullfile(xsd_path, 'SICD_schema*.xsd'));
for i = 1:numel(filelist) % Try with all known XSDs
    schema = factory.newSchema( java.io.File(fullfile(xsd_path, filelist(i).name)) );
    try
        schema.newValidator().validate( javax.xml.transform.stream.StreamSource( java.io.StringBufferInputStream( SICD_xml_string ) ) );
        disp(['XML passed validation for ' filelist(i).name '.']);
        validation_report = cell(0,3); % Clear any error messages from earlier schemas not validating
        if i ~= numel(filelist) % Validation is not against most recent SICD version
            % Report if schema validated against is not the most current.
            validation_report = add_val_inc(validation_report, 'Warning', ...
                'SICD version not current.', ...
                ['SICD schema validated: ' filelist(i).name, '.  Most recent schema: ' filelist(end).name]);
        end
        break;
    catch ME
        % Find message after "Java exception occurred", but before stack trace
        ind = find(ME.message<=13,2,'first'); % Probably a cleaner way with regexp
        validation_report = add_val_inc(validation_report, 'Error', ...
            ['XML failed validation for ' filelist(i).name '.'], ...
            ['First reason: ' ME.message((ind(1)+1):(ind(2)-1))]);
    end
end

try % If schema validation fails, no guarantee that any of the following will work
%% 2. Next check for internal consistency of data
% 2.1. Scalar TimeCOAPoly means SPOTLIGHT data
if strcmp(SICD_meta.CollectionInfo.RadarMode.ModeType,'SPOTLIGHT')
    if any(SICD_meta.Grid.TimeCOAPoly(2:end)~=0)
        validation_report = add_val_inc(validation_report, 'Error', ...
            'SPOTLIGHT data should only have scalar TimeCOAPoly.', ...
            ['TimeCOAPoly: ' mat2str(SICD_meta.Grid.TimeCOAPoly)]);
    end
else % Non-spotlight data
    if all(SICD_meta.Grid.TimeCOAPoly(2:end)==0)
        validation_report = add_val_inc(validation_report, 'Warning', ...
            'Non-SPOTLIGHT data will generally have more than one nonzero term in TimeCOAPoly unless "formed as spotlight".', ...
            ['Grid.TimeCOAPoly: ' mat2str(SICD_meta.Grid.TimeCOAPoly)]);
    end
end
% 2.2. FFT signs in both dimensions almost certainly have to be equal
if SICD_meta.Grid.Col.Sgn~=SICD_meta.Grid.Row.Sgn
    validation_report = add_val_inc(validation_report, 'Error', ...
        'FFT signs in row and column direction should be the same.', ...
        ['Grid.Row.Sgn: ' num2str(SICD_meta.Grid.Row.Sgn) ...
        ' , Grid.Col.Sgn: ' num2str(SICD_meta.Grid.Col.Sgn)]);
end

% 2.3. Frequency support parameters
% We add the eps in all these tests to handle numerical precision
% issues for two values that are very nearly equal.
SF_BOUNDS_STR = 'Violation of spatial frequency extent bounds.';
% 2.3.1.
if SICD_meta.Grid.Col.DeltaK2<=SICD_meta.Grid.Col.DeltaK1
    validation_report = add_val_inc(validation_report, 'Error', ...
        SF_BOUNDS_STR, ...
        ['Grid.Col.DeltaK1: ' num2str(SICD_meta.Grid.Col.DeltaK1) ...
        ' , Grid.Col.DeltaK2: ' num2str(SICD_meta.Grid.Col.DeltaK2)]);
else
    % 2.3.2.
    if SICD_meta.Grid.Col.DeltaK2>((1/(2*SICD_meta.Grid.Col.SS))+eps(SICD_meta.Grid.Col.SS))
        validation_report = add_val_inc(validation_report, 'Error', ...
            SF_BOUNDS_STR, ...
            ['0.5/Grid.Col.SS: ' num2str(0.5/SICD_meta.Grid.Col.SS) ...
            ' , Grid.Col.DeltaK2: ' num2str(SICD_meta.Grid.Col.DeltaK2)]);
    end
    % 2.3.3.
    if SICD_meta.Grid.Col.DeltaK1<((-1/(2*SICD_meta.Grid.Col.SS))-eps(SICD_meta.Grid.Col.SS))
        validation_report = add_val_inc(validation_report, 'Error', ...
            SF_BOUNDS_STR, ...
            ['0.5/Grid.Col.SS: ' num2str(0.5/SICD_meta.Grid.Col.SS) ...
            ' , Grid.Col.DeltaK1: ' num2str(SICD_meta.Grid.Col.DeltaK1)]);
    end
    % 2.3.4.
    if SICD_meta.Grid.Col.ImpRespBW>((SICD_meta.Grid.Col.DeltaK2-...
            SICD_meta.Grid.Col.DeltaK1)+eps(SICD_meta.Grid.Col.ImpRespBW))
        validation_report = add_val_inc(validation_report, 'Error', ...
            SF_BOUNDS_STR, ...
            ['Grid.Col.ImpRespBW: ' num2str(SICD_meta.Grid.Col.ImpRespBW) ...
            ' , Grid.Col.DeltaK2-Grid.Col.DeltaK1: ' ...
            num2str(SICD_meta.Grid.Col.DeltaK2-SICD_meta.Grid.Col.DeltaK1)]);
    end
end
% 2.3.5.
if SICD_meta.Grid.Row.DeltaK2<=SICD_meta.Grid.Row.DeltaK1
    validation_report = add_val_inc(validation_report, 'Error', ...
        SF_BOUNDS_STR, ...
        ['Grid.Row.DeltaK1: ' num2str(SICD_meta.Grid.Row.DeltaK1) ...
        ' , Grid.Row.DeltaK2: ' num2str(SICD_meta.Grid.Row.DeltaK2)]);
else
    % 2.3.6.
    if SICD_meta.Grid.Row.DeltaK2>((1/(2*SICD_meta.Grid.Row.SS))+eps(SICD_meta.Grid.Row.SS))
        validation_report = add_val_inc(validation_report, 'Error', ...
            SF_BOUNDS_STR, ...
            ['0.5/Grid.Row.SS: ' num2str(0.5/SICD_meta.Grid.Row.SS) ...
            ' , Grid.Col.DeltaK2: ' num2str(SICD_meta.Grid.Row.DeltaK2)]);
    end
    % 2.3.7.
    if SICD_meta.Grid.Row.DeltaK1<((-1/(2*SICD_meta.Grid.Row.SS))-eps(SICD_meta.Grid.Row.SS))
        validation_report = add_val_inc(validation_report, 'Error', ...
            SF_BOUNDS_STR, ...
            ['0.5/Grid.Row.SS: ' num2str(0.5/SICD_meta.Grid.Row.SS) ...
            ' , Grid.Row.DeltaK1: ' num2str(SICD_meta.Grid.Row.DeltaK1)]);
    end
    % 2.3.8.
    if SICD_meta.Grid.Row.ImpRespBW>((SICD_meta.Grid.Row.DeltaK2-...
            SICD_meta.Grid.Row.DeltaK1)+eps(SICD_meta.Grid.Row.ImpRespBW))
        validation_report = add_val_inc(validation_report, 'Error', ...
            SF_BOUNDS_STR, ...
            ['Grid.Row.ImpRespBW: ' num2str(SICD_meta.Grid.Row.ImpRespBW) ...
            ' , Grid.Row.DeltaK2-Grid.Row.DeltaK1: ' ...
            num2str(SICD_meta.Grid.Row.DeltaK2-SICD_meta.Grid.Row.DeltaK1)]);
    end
    % 2.3.9. Compute our own DeltaK1/K2 and test for consistency with
    % DeltaKCOAPoly, ImpRespBW, and SS.  Here, we assume the min and max of
    % DeltaKCOAPoly must be on the vertices of the image, since it is
    % smooth and monotonic in most cases-- although in actuality this is
    % not always the case.  To be totally generic, we would have to search
    % for an interior min and max as well.
    if isfield(SICD_meta.ImageData,'ValidData') % Test vertices
        vertices = double([SICD_meta.ImageData.ValidData.Vertex.Col;...
            SICD_meta.ImageData.ValidData.Vertex.Row]);
    else % Use edges of full image
        vertices = [0 SICD_meta.ImageData.NumCols-1 SICD_meta.ImageData.NumCols-1 0; ...
            0 0 SICD_meta.ImageData.NumRows-1 SICD_meta.ImageData.NumRows-1];
    end
    vertices = single(vertices); % SICD stores indices as ints
    if isfield(SICD_meta.Grid.Row,'DeltaKCOAPoly')
        min_row_dk = Inf; max_row_dk= -Inf;
    else
        min_row_dk = 0; max_row_dk= 0;        
    end
    if isfield(SICD_meta.Grid.Col,'DeltaKCOAPoly')
        min_col_dk = Inf; max_col_dk = -Inf;
    else
        min_col_dk = 0; max_col_dk = 0;
    end
    for i = 1:size(vertices,2)
        if isfield(SICD_meta.Grid.Row,'DeltaKCOAPoly')
            currentDeltaK = sicd_polyval2d(SICD_meta.Grid.Row.DeltaKCOAPoly,...
                vertices(1,i),vertices(2,i),SICD_meta);
            min_row_dk = min(min_row_dk, currentDeltaK);
            max_row_dk = max(max_row_dk, currentDeltaK);
        end
        if isfield(SICD_meta.Grid.Col,'DeltaKCOAPoly')
            currentDeltaK = sicd_polyval2d(SICD_meta.Grid.Col.DeltaKCOAPoly,...
                vertices(1,i),vertices(2,i),SICD_meta);
            min_col_dk = min(min_col_dk, currentDeltaK);
            max_col_dk = max(max_col_dk, currentDeltaK);
        end
    end
    min_row_dk = min_row_dk - (SICD_meta.Grid.Row.ImpRespBW/2);
    max_row_dk = max_row_dk + (SICD_meta.Grid.Row.ImpRespBW/2);
    min_col_dk = min_col_dk - (SICD_meta.Grid.Col.ImpRespBW/2);
    max_col_dk = max_col_dk + (SICD_meta.Grid.Col.ImpRespBW/2);
    % Wrapped spectrum
    if (min_row_dk < -(1/SICD_meta.Grid.Row.SS)/2) || ...
            (max_row_dk > (1/SICD_meta.Grid.Row.SS)/2)
        min_row_dk = -(1/SICD_meta.Grid.Row.SS)/2;
        max_row_dk = -min_row_dk;
    end
    if (min_col_dk < -(1/SICD_meta.Grid.Col.SS)/2) || ...
            (max_col_dk > (1/SICD_meta.Grid.Col.SS)/2)
        min_col_dk = -(1/SICD_meta.Grid.Col.SS)/2;
        max_col_dk = -min_col_dk;
    end
    DK_TOL = 1e-2;
    % 2.3.9.1
    if abs(SICD_meta.Grid.Row.DeltaK1/min_row_dk - 1) > DK_TOL
        validation_report = add_val_inc(validation_report, 'Error', ...
            SF_BOUNDS_STR, ...
            ['Grid.Row.DeltaK1: ' num2str(SICD_meta.Grid.Row.DeltaK1) ...
            ' , Derived Grid.Row.DeltaK1: ' num2str(min_row_dk)]);
    end
    % 2.3.9.2
    if abs(SICD_meta.Grid.Row.DeltaK2/max_row_dk - 1) > DK_TOL
        validation_report = add_val_inc(validation_report, 'Error', ...
            SF_BOUNDS_STR, ...
            ['Grid.Row.DeltaK2: ' num2str(SICD_meta.Grid.Row.DeltaK2) ...
            ' , Derived Grid.Row.DeltaK2: ' num2str(max_row_dk)]);
    end
    % 2.3.9.3
    if abs(SICD_meta.Grid.Col.DeltaK1/min_col_dk - 1) > DK_TOL
        validation_report = add_val_inc(validation_report, 'Error', ...
            SF_BOUNDS_STR, ...
            ['Grid.Col.DeltaK1: ' num2str(SICD_meta.Grid.Col.DeltaK1) ...
            ' , Derived Grid.Col.DeltaK1: ' num2str(min_col_dk)]);
    end
    % 2.3.9.4
    if abs(SICD_meta.Grid.Col.DeltaK2/max_col_dk - 1) > DK_TOL
        validation_report = add_val_inc(validation_report, 'Error', ...
            SF_BOUNDS_STR, ...
            ['Grid.Col.DeltaK2: ' num2str(SICD_meta.Grid.Col.DeltaK2) ...
            ' , Derived Grid.Col.DeltaK2: ' num2str(max_col_dk)]);
    end
end
if isfield(SICD_meta, 'PFA') % A further check for PFA data
    % Slow-time deskew would allow for PFA.Kaz2-PFA.Kaz1>(1/Grid.Col.SS),
    % since Kaz bandwidth is compressed from original polar annulus.
    if ~isfield(SICD_meta.PFA,'STDeskew') || ...
            ~SICD_meta.PFA.STDeskew.Applied 
        % 2.3.10.
        if (SICD_meta.PFA.Kaz2-SICD_meta.Grid.Col.KCtr)>...
                ((1/(2*SICD_meta.Grid.Col.SS))+eps(SICD_meta.Grid.Col.SS))
            validation_report = add_val_inc(validation_report, 'Error', ...
                SF_BOUNDS_STR, ...
                ['0.5/Grid.Col.SS: ' num2str(0.5/SICD_meta.Grid.Col.SS) ...
                ' , PFA.Kaz2-Grid.Col.KCtr: ' ...
                num2str(SICD_meta.PFA.Kaz2-SICD_meta.Grid.Col.KCtr)]);
        end
        % 2.3.11.
        if (SICD_meta.PFA.Kaz1-SICD_meta.Grid.Col.KCtr)<...
                ((-1/(2*SICD_meta.Grid.Col.SS))-eps(SICD_meta.Grid.Col.SS))
            validation_report = add_val_inc(validation_report, 'Error', ...
                SF_BOUNDS_STR, ...
                ['0.5/Grid.Col.SS: ' num2str(0.5/SICD_meta.Grid.Col.SS) ...
                ' , PFA.Kaz1-Grid.Col.KCtr: ' ...
                num2str(SICD_meta.PFA.Kaz1-SICD_meta.Grid.Col.KCtr)]);
        end
    end
    % 2.3.12.
    if (SICD_meta.PFA.Krg2-SICD_meta.Grid.Row.KCtr)>...
            ((1/(2*SICD_meta.Grid.Row.SS))+eps(SICD_meta.Grid.Row.SS))
        validation_report = add_val_inc(validation_report, 'Error', ...
            SF_BOUNDS_STR, ...
            ['0.5/Grid.Row.SS: ' num2str(0.5/SICD_meta.Grid.Row.SS) ...
            ' , PFA.Krg2-Grid.Row.KCtr: ' ...
            num2str(SICD_meta.PFA.Krg2-SICD_meta.Grid.Row.KCtr)]);
    end
    % 2.3.13.
    if (SICD_meta.PFA.Krg1-SICD_meta.Grid.Row.KCtr)<...
            ((-1/(2*SICD_meta.Grid.Row.SS))-eps(SICD_meta.Grid.Row.SS))
        validation_report = add_val_inc(validation_report, 'Error', ...
            SF_BOUNDS_STR, ...
            ['0.5/Grid.Row.SS: ' num2str(0.5/SICD_meta.Grid.Row.SS) ...
            ' , PFA.Krg1-Grid.Row.KCtr: ' ...
            num2str(SICD_meta.PFA.Krg1-SICD_meta.Grid.Row.KCtr)]);
    end
    % 2.3.14.
    if SICD_meta.Grid.Col.ImpRespBW>((SICD_meta.PFA.Kaz2-...
            SICD_meta.PFA.Kaz1)+eps(SICD_meta.Grid.Col.ImpRespBW))
        validation_report = add_val_inc(validation_report, 'Error', ...
            SF_BOUNDS_STR, ...
            ['Grid.Col.ImpRespBW: ' num2str(SICD_meta.Grid.Col.ImpRespBW) ...
            ' , PFA.Kaz2-PFA.Kaz1: ' num2str(SICD_meta.PFA.Kaz2-SICD_meta.PFA.Kaz1)]);
    end
    % 2.3.15.
    if SICD_meta.Grid.Row.ImpRespBW>((SICD_meta.PFA.Krg2-...
            SICD_meta.PFA.Krg1)+eps(SICD_meta.Grid.Row.ImpRespBW))
        validation_report = add_val_inc(validation_report, 'Error', ...
            SF_BOUNDS_STR, ...
            ['Grid.Row.ImpRespBW: ' num2str(SICD_meta.Grid.Row.ImpRespBW) ...
            ' , PFA.Krg2-PFA.Krg1: ' num2str(SICD_meta.PFA.Krg2-SICD_meta.PFA.Krg1)]);
    end
    % 2.3.16.
    if (SICD_meta.Grid.Col.KCtr ~= 0) && ...
       (abs(SICD_meta.Grid.Col.KCtr - mean([SICD_meta.PFA.Kaz1 SICD_meta.PFA.Kaz2])) > 1e-5)
        validation_report = add_val_inc(validation_report, 'Warning', ...
            SF_BOUNDS_STR, ...
            ['Grid.Col.KCtr: ' num2str(SICD_meta.Grid.Col.KCtr) ...
            ' , mean([PFA.Kaz2, PFA.Kaz1]): ' ...
            num2str(mean([SICD_meta.PFA.Kaz1 SICD_meta.PFA.Kaz2]))]);
    end
end

% 2.4. Does WgtFunct agree with WgtType?
% Sloppy code, since multiple lines are cut-and-paste for row and column,
% rather than wrapped as a function.
DEFAULT_WGT_SIZE = 512;
WGT_TOL = 1e-3;
sicd_meta_temp = SICD_meta; % We may need to take out WgtFunct to repopulate it
% 2.4.1
if isfield(SICD_meta.Grid.Row,'WgtFunct') % WgtFunct is optional field
    sicd_meta_temp.Grid.Row = rmfield(sicd_meta_temp.Grid.Row,'WgtFunct');
    rowWgtFunct = SICD_meta.Grid.Row.WgtFunct; % Use WgtFunct if it exist for further functions
end
if isfield(SICD_meta.Grid.Row,'WgtType') % WgtType is optional field
    try
        rowFun = sicdweight2fun(sicd_meta_temp.Grid.Row); % Will error if WgtFunct cannot be created from WgtType
        if isfield(SICD_meta.Grid.Row,'WgtFunct') % Optional field
            numwgts = numel(SICD_meta.Grid.Row.WgtFunct);
            if (isempty(rowFun) && ...% sicdweight2fun way of passing UNIFORM
                    any(diff(SICD_meta.Grid.Row.WgtFunct))) || ...
                    (max(abs(SICD_meta.Grid.Row.WgtFunct - ... % Non-uniform case
                    rowFun(numwgts))) > WGT_TOL)
                validation_report = add_val_inc(validation_report, 'Error', ...
                    'Row.WgtFunct values inconsistent with Row.WgtType', ...
                    ['Grid.Row.WgtType.WindowName: ' SICD_meta.Grid.Row.WgtType.WindowName]);
            end
        elseif ~isempty(rowFun)
            rowWgtFunct = rowFun(DEFAULT_WGT_SIZE); % If no WgtFunct, use WgtType computation for future tests
        else
            rowWgtFunct = []; % UNIFORM weighting is special case
        end
    catch
        validation_report = add_val_inc(validation_report, 'Style', ...
            'Unrecognized weighting description', ...
            ['Grid.Row.WgtType.WindowName: ' SICD_meta.Grid.Row.WgtType.WindowName]);
    end
end
% 2.4.2
if isfield(SICD_meta.Grid.Col,'WgtFunct') % WgtFunct is optional field
    sicd_meta_temp.Grid.Col = rmfield(sicd_meta_temp.Grid.Col,'WgtFunct');
    colWgtFunct = SICD_meta.Grid.Col.WgtFunct; % Use WgtFunct if it exist for further functions
end
if isfield(SICD_meta.Grid.Col,'WgtType') % WgtType is optional field
    try
        colFun = sicdweight2fun(sicd_meta_temp.Grid.Col); % Will error if WgtFunct cannot be created from WgtType
        if isfield(SICD_meta.Grid.Col,'WgtFunct') % Optional field
            numwgts = numel(SICD_meta.Grid.Col.WgtFunct);
            if (isempty(colFun) && ...% sicdweight2fun way of passing UNIFORM
                    any(diff(SICD_meta.Grid.Col.WgtFunct))) || ...
                    (max(abs(SICD_meta.Grid.Col.WgtFunct - ... % Non-uniform case
                    colFun(numwgts))) > WGT_TOL)
                validation_report = add_val_inc(validation_report, 'Error', ...
                    'Col.WgtFunct values inconsistent with Col.WgtType', ...
                    ['Grid.Col.WgtType.WindowName: ' SICD_meta.Grid.Col.WgtType.WindowName]);
            end
        elseif ~isempty(colFun)
            colWgtFunct = colFun(DEFAULT_WGT_SIZE); % If no WgtFunct, use WgtType computation for future tests
        else
            colWgtFunct = []; % UNIFORM weighting is special case
        end
    catch
        validation_report = add_val_inc(validation_report, 'Style', ...
            'Unrecognized weighting description', ...
            ['Grid.Col.WgtType.WindowName: ' SICD_meta.Grid.Col.WgtType.WindowName]);
    end
end
% 2.4.3
if isfield(SICD_meta.Grid.Col,'WgtType') && ...
        isfield(SICD_meta.Grid.Col.WgtType,'WindowName') && ...
        ~any(strcmp(SICD_meta.Grid.Col.WgtType.WindowName,{'UNIFORM','UNKNOWN'})) && ...
        ~isfield(SICD_meta.Grid.Col,'WgtFunct')
    validation_report = add_val_inc(validation_report, 'Warning', ...
        'Non-uninform weighting, but no WgtFunct provided.', ...
        ['Grid.Col.WgtType.WindowName: ' SICD_meta.Grid.Col.WgtType.WindowName ]);
end
% 2.4.4
if isfield(SICD_meta.Grid.Row,'WgtType') && ...
        isfield(SICD_meta.Grid.Row.WgtType,'WindowName') && ...
        ~any(strcmp(SICD_meta.Grid.Row.WgtType.WindowName,{'UNIFORM','UNKNOWN'})) && ...
        ~isfield(SICD_meta.Grid.Row,'WgtFunct')
    validation_report = add_val_inc(validation_report, 'Warning', ...
        'Non-uninform weighting, but no WgtFunct provided.', ...
        ['Grid.Row.WgtType.WindowName: ' SICD_meta.Grid.Row.WgtType.WindowName ]);
end
% 2.5. Is ImpRespWid consistent with WgtFunct/WgtType?
IPR_TOL = 1e-2; % ImpRespWid should agree within this tolerance (relative to computed IPR)
if exist('rowWgtFunct','var') % From optional fields
    der_ImpRespWid = computeImpRespWid(rowWgtFunct, SICD_meta.Grid.Row.ImpRespBW);
    if abs(SICD_meta.Grid.Row.ImpRespWid - der_ImpRespWid) > (IPR_TOL*der_ImpRespWid)
        validation_report = add_val_inc(validation_report, 'Error', ...
            'Grid.Row.ImpRespWid not consistent with other Grid.Row data', ...
            ['Grid.Row.ImpRespWid: ' num2str(SICD_meta.Grid.Row.ImpRespWid) ...
            ' , Derived Row.ImpRespWid: ' num2str(der_ImpRespWid) ]);
    end
end
if exist('colWgtFunct','var') % From optional fields
    der_ImpRespWid = computeImpRespWid(colWgtFunct, SICD_meta.Grid.Col.ImpRespBW);
    if abs(SICD_meta.Grid.Col.ImpRespWid - der_ImpRespWid) > (IPR_TOL*der_ImpRespWid)
        validation_report = add_val_inc(validation_report, 'Error', ...
            'Grid.Col.ImpRespWid not consistent with other Grid.Col data', ...
            ['Grid.Col.ImpRespWid: ' num2str(SICD_meta.Grid.Col.ImpRespWid) ...
            ' , Derived Col.ImpRespWid: ' num2str(der_ImpRespWid) ]);
    end
end

% 2.6
ARPPoly_order = max([find(SICD_meta.Position.ARPPoly.X,1,'last'), ...
    find(SICD_meta.Position.ARPPoly.Y,1,'last'), ...
    find(SICD_meta.Position.ARPPoly.Z,1,'last')]);
if ARPPoly_order<2 % Must be able to derive at least positon and velocity
    validation_report = add_val_inc(validation_report, 'Error', ...
        'ARPPoly should have at least position and velocity terms.', ...
        ['Position.ARPPoly.X: ' mat2str(SICD_meta.Position.ARPPoly.X,3) ...
        ' , Position.ARPPoly.Y: ' mat2str(SICD_meta.Position.ARPPoly.Y,3) ...
        ' , Position.ARPPoly.Z: ' mat2str(SICD_meta.Position.ARPPoly.Z,3)]);
end
% 2.7 Assure derived fields in SCPCOA are consistent with other data
sicd_meta_temp = SICD_meta; % We may need to take out WgtFunct to repopulate it
sicd_meta_temp = rmfield(sicd_meta_temp,'SCPCOA');
sicd_meta_temp = derived_sicd_fields(sicd_meta_temp); % Re-derive all SCPCOA values
SCPCOA_INCONSISTENT_STR = 'SCPCOA fields not consistent with other SICD metadata.';
SCPCOA_TOL = 1e-2;
% 2.7.1 Special test for SCPTime with stricter tolerance because it should
% be exact and everything else depends on it.
if abs(SICD_meta.SCPCOA.SCPTime - ...
        SICD_meta.Grid.TimeCOAPoly(1)) > eps(SICD_meta.Grid.TimeCOAPoly(1))
    validation_report = add_val_inc(validation_report, 'Error', ...
        SCPCOA_INCONSISTENT_STR, ...
        ['SCPCOA.SCPTime - Grid.TimeCOAPoly(1): ' ...
        num2str(SICD_meta.SCPCOA.SCPTime - SICD_meta.Grid.TimeCOAPoly(1)) ]);
end
% 2.7.2-2.7.4 XYZ fields of SCPCOA
SCPCOA_XYZ_FIELDS = {'ARPPos', 'ARPVel'}; % Could check 'ARPAcc' as well,
% but ARPAcc doesn't really matter and the instantaneous ARPAcc can often
% vary from the ARPPoly polynomial derived one.
for i = 1:numel(SCPCOA_XYZ_FIELDS)
    pos_diff = norm([SICD_meta.SCPCOA.(SCPCOA_XYZ_FIELDS{i}).X - ...
        sicd_meta_temp.SCPCOA.(SCPCOA_XYZ_FIELDS{i}).X ...
        SICD_meta.SCPCOA.(SCPCOA_XYZ_FIELDS{i}).Y - ...
        sicd_meta_temp.SCPCOA.(SCPCOA_XYZ_FIELDS{i}).Y ...
        SICD_meta.SCPCOA.(SCPCOA_XYZ_FIELDS{i}).Z - ...
        sicd_meta_temp.SCPCOA.(SCPCOA_XYZ_FIELDS{i}).Z]);
    if pos_diff > SCPCOA_TOL
        validation_report = add_val_inc(validation_report, 'Error', ...
            SCPCOA_INCONSISTENT_STR, ...
            ['SCPCOA.' (SCPCOA_XYZ_FIELDS{i}) ' - Derived SCPCOA.' ...
            (SCPCOA_XYZ_FIELDS{i}) ': ' num2str(pos_diff)]);
    end
end
% 2.7.5 SideOfTrack
if ~strcmp(SICD_meta.SCPCOA.SideOfTrack, sicd_meta_temp.SCPCOA.SideOfTrack)
    validation_report = add_val_inc(validation_report, 'Error', ...
        SCPCOA_INCONSISTENT_STR, ...
        ['SCPCOA.SideOfTrack: ' SICD_meta.SCPCOA.SideOfTrack ...
        ' , Derived SCPCOA.SideOfTrack: ' sicd_meta_temp.SCPCOA.SideOfTrack ]);
end
% 2.7.6-2.7.14 Scalar numeric fields of SCPCOA
SCPCOA_NUM_FIELDS = {'SlantRange', 'GroundRange', ...
    'DopplerConeAng', 'GrazeAng', 'IncidenceAng', 'TwistAng', ...
    'SlopeAng', 'AzimAng', 'LayoverAng'};
for i = 1:numel(SCPCOA_NUM_FIELDS)
    if abs(SICD_meta.SCPCOA.(SCPCOA_NUM_FIELDS{i}) - ...
            sicd_meta_temp.SCPCOA.(SCPCOA_NUM_FIELDS{i})) > SCPCOA_TOL
        validation_report = add_val_inc(validation_report, 'Error', ...
            SCPCOA_INCONSISTENT_STR, ...
            ['SCPCOA.' SCPCOA_NUM_FIELDS{i} ': ' ...
            num2str(SICD_meta.SCPCOA.(SCPCOA_NUM_FIELDS{i})) ...
            ' , Derived SCPCOA.' SCPCOA_NUM_FIELDS{i} ': ' ...
            num2str(sicd_meta_temp.SCPCOA.(SCPCOA_NUM_FIELDS{i})) ]);
    end
end

% 2.8 Waveform description consistency
WF_TOL = 1e-3; % Relative tolerance
WF_INCONSISTENT_STR = 'Waveform fields not consistent.';
try % All WFParameters fields are optional
    wfmin = min([SICD_meta.RadarCollection.Waveform.WFParameters.TxFreqStart]);
    wfmax = max([SICD_meta.RadarCollection.Waveform.WFParameters.TxFreqStart] + ...
        [SICD_meta.RadarCollection.Waveform.WFParameters.TxRFBandwidth]);
end
% 2.8.1
try % All WFParameters fields are optional
    if abs((wfmin/SICD_meta.RadarCollection.TxFrequency.Min) - 1) > WF_TOL
        validation_report = add_val_inc(validation_report, 'Error', ...
            WF_INCONSISTENT_STR, ...
            ['RadarCollection.Waveform.WFParameters.TxFreqStart: ' ...
            num2str(wfmin) ...
            ' , RadarCollection.TxFrequency.Min: ' ...
            num2str(SICD_meta.RadarCollection.TxFrequency.Min) ]);
    end
end
% 2.8.2
try % All WFParameters fields are optional
    if abs((wfmax/SICD_meta.RadarCollection.TxFrequency.Max) - 1) > WF_TOL
        validation_report = add_val_inc(validation_report, 'Error', ...
            WF_INCONSISTENT_STR, ...
            ['RadarCollection.Waveform.WFParameters.TxFreqStart+TxRFBandwidth: ' ...
            num2str(wfmax) ...
            ' , RadarCollection.TxFrequency.Max: ' ...
            num2str(SICD_meta.RadarCollection.TxFrequency.Max)]);
    end
end
if isfield(SICD_meta.RadarCollection,'Waveform')
    for i = 1:numel(SICD_meta.RadarCollection.Waveform.WFParameters)
        wfp = SICD_meta.RadarCollection.Waveform.WFParameters(i); % To shorten notation
        % 2.8.3
        try % All WFParameters fields are optional
            if abs((wfp.TxRFBandwidth/(wfp.TxPulseLength * wfp.TxFMRate)) - 1) > WF_TOL
                validation_report = add_val_inc(validation_report, 'Error', ...
                    WF_INCONSISTENT_STR, ...
                    ['RadarCollection.Waveform.WFParameters.TxRFBandwidth: ' ...
                    num2str(wfp.TxRFBandwidth) ...
                    ' , RadarCollection.TxFrequency.TxFMRate*TxPulseLength: ' ...
                    num2str(wfp.TxFMRate * ...
                    wfp.TxPulseLength) ]);
            end
        end
        % 2.8.4
        try % All WFParameters fields are optional
            if strcmp(wfp.RcvDemodType, 'CHIRP') && (wfp.RcvFMRate ~= 0)
                validation_report = add_val_inc(validation_report, 'Error', ...
                    WF_INCONSISTENT_STR, ...
                    ['RadarCollection.Waveform.WFParameters.RcvDemodType: ' ...
                    wfp.RcvDemodType ...
                    ' , RadarCollection.Waveform.WFParameters.RcvFMRate: ' ...
                    num2str(wfp.RcvFMRate) ]);
            end
        end
        % 2.8.5
        try % All WFParameters fields are optional
            if strcmp(wfp.RcvDemodType, 'STRETCH') && ...
                    abs((wfp.RcvFMRate/wfp.TxFMRate) - 1) > WF_TOL
                validation_report = add_val_inc(validation_report, 'Warning', ... % Physically possible but unlikely
                    WF_INCONSISTENT_STR, ...
                    ['RadarCollection.Waveform.WFParameters.RcvDemodType: ' ...
                    wfp.RcvDemodType ...
                    ' , RadarCollection.Waveform.WFParameters.RcvFMRate: ' ...
                    num2str(wfp.RcvFMRate) ...
                    ' , RadarCollection.Waveform.WFParameters.TxFMRate: ' ...
                    num2str(wfp.TxFMRate)]);
            end
        end
        % 2.8.6
        % Absolute frequencies must be positive
        if ~isfield(SICD_meta.RadarCollection,'RefFreqIndex') && ...
                SICD_meta.RadarCollection.TxFrequency.Min<=0
            validation_report = add_val_inc(validation_report, 'Error', ...
                WF_INCONSISTENT_STR, ...
                ['RadarCollection.TxFrequency.Min: ' ...
                SICD_meta.RadarCollection.TxFrequency.Min]);
        end
        % 2.8.7
        try % All WFParameters fields are optional
            % Absolute frequencies must be positive
            if ~isfield(SICD_meta.RadarCollection,'RefFreqIndex') && ...
                    wfp.TxFreqStart<=0
                validation_report = add_val_inc(validation_report, 'Error', ...
                    WF_INCONSISTENT_STR, ...
                    ['RadarCollection.Waveform.WFParameters.TxFreqStart: ' ...
                    wfp.TxFreqStart]);
            end
        end
        % 2.8.8
        try % All WFParameters fields are optional
            % Absolute frequencies must be positive
            if ~isfield(SICD_meta.RadarCollection,'RefFreqIndex') && ...
                    wfp.RcvFreqStart<=0
                validation_report = add_val_inc(validation_report, 'Error', ...
                    WF_INCONSISTENT_STR, ...
                    ['RadarCollection.Waveform.WFParameters.RcvFreqStart: ' ...
                    wfp.RcvFreqStart]);
            end
        end
        % 2.8.9
        try % All WFParameters fields are optional
            if wfp.TxPulseLength > wfp.RcvWindowLength
                validation_report = add_val_inc(validation_report, 'Error', ...
                    WF_INCONSISTENT_STR, ...
                    ['RadarCollection.Waveform.WFParameters.TxPulseLength: ' ...
                    num2str(wfp.TxPulseLength) ...
                    ' , RadarCollection.Waveform.WFParameters.RcvWindowLength: ' ...
                    num2str(wfp.RcvWindowLength) ]);
            end
        end
        % 2.8.10
        try % All WFParameters fields are optional
            if wfp.RcvIFBandwidth > wfp.ADCSampleRate
                validation_report = add_val_inc(validation_report, 'Error', ...
                    WF_INCONSISTENT_STR, ...
                    ['RadarCollection.Waveform.WFParameters.RcvIFBandwidth: ' ...
                    num2str(wfp.RcvIFBandwidth) ...
                    ' , RadarCollection.Waveform.WFParameters.ADCSampleRate: ' ...
                    num2str(SICD_meta.RadarCollection.Waveform.WFParameters.ADCSampleRate) ]);
            end
        end
        % 2.8.11
        try % All WFParameters fields are optional
            if strcmp(wfp.RcvDemodType, 'CHIRP') && ...
                    (wfp.TxRFBandwidth > wfp.ADCSampleRate)
                validation_report = add_val_inc(validation_report, 'Error', ...
                    WF_INCONSISTENT_STR, ...
                    ['RadarCollection.Waveform.WFParameters.RcvDemodType: ' ...
                    wfp.RcvDemodType ...
                    ' , RadarCollection.Waveform.WFParameters.TxRFBandwidth: ' ...
                    num2str(wfp.TxRFBandwidth) ...
                    ' , RadarCollection.Waveform.WFParameters.ADCSampleRate: ' ...
                    num2str(wfp.ADCSampleRate)]);
            end
        end
        % 2.8.12
        try % All WFParameters fields are optional
            freq_tol = (wfp.RcvWindowLength-wfp.TxPulseLength) * wfp.TxFMRate;
            if (wfp.RcvFreqStart >= (wfp.TxFreqStart+wfp.TxRFBandwidth+freq_tol)) || ...
                    (wfp.RcvFreqStart <= wfp.TxFreqStart-freq_tol)
                validation_report = add_val_inc(validation_report, 'Error', ...
                    WF_INCONSISTENT_STR, ...
                    ['RadarCollection.Waveform.WFParameters.RcvFreqStart: ' ...
                    num2str(wfp.RcvFreqStart) ]);
            end
        end
    end
end
% 2.8.13
if ~(SICD_meta.RadarCollection.TxFrequency.Min<...
        SICD_meta.RadarCollection.TxFrequency.Max)
    validation_report = add_val_inc(validation_report, 'Error', ...
        'Transmit frequency bounds not valid.', ...
        ['SICD_meta.RadarCollection.TxFrequency.Min = ' ...
        num2str(SICD_meta.RadarCollection.TxFrequency.Min) ...
        ', SICD_meta.RadarCollection.TxFrequency.Max: ' ...
        num2str(SICD_meta.RadarCollection.TxFrequency.Max)]);
end
% 2.8.14
if ~(SICD_meta.ImageFormation.TxFrequencyProc.MinProc<...
        SICD_meta.ImageFormation.TxFrequencyProc.MaxProc)
    validation_report = add_val_inc(validation_report, 'Error', ...
        'Image formation frequency bounds not valid.', ...
        ['SICD_meta.ImageFormation.TxFrequencyProc.MinProc = ' ...
        num2str(SICD_meta.ImageFormation.TxFrequencyProc.MinProc) ...
        ', SICD_meta.ImageFormation.TxFrequencyProc.MaxProc: ' ...
        num2str(SICD_meta.ImageFormation.TxFrequencyProc.MaxProc)]);
end

% 2.9 Image corners coordinate ordering
% The corner coordinates are allowed to be an approximation by the spec, so
% we don't check for precision here (one also doesn't know the HAE used to
% compute Lat/Lon for each corner), but we do at least check that the
% corners are ordered properly.
corner_latlons = [SICD_meta.GeoData.ImageCorners.ICP.FRFC.Lat ...
    SICD_meta.GeoData.ImageCorners.ICP.FRLC.Lat ...
    SICD_meta.GeoData.ImageCorners.ICP.LRLC.Lat ...
    SICD_meta.GeoData.ImageCorners.ICP.LRFC.Lat; ...
    SICD_meta.GeoData.ImageCorners.ICP.FRFC.Lon ...
    SICD_meta.GeoData.ImageCorners.ICP.FRLC.Lon ...
    SICD_meta.GeoData.ImageCorners.ICP.LRLC.Lon ...
    SICD_meta.GeoData.ImageCorners.ICP.LRFC.Lon];
new_corners_meta = add_sicd_corners(SICD_meta); % Recompute corner coordinates from SICD sensor model
new_corner_latlons = [new_corners_meta.GeoData.ImageCorners.ICP.FRFC.Lat ...
    new_corners_meta.GeoData.ImageCorners.ICP.FRLC.Lat ...
    new_corners_meta.GeoData.ImageCorners.ICP.LRLC.Lat ...
    new_corners_meta.GeoData.ImageCorners.ICP.LRFC.Lat; ...
    new_corners_meta.GeoData.ImageCorners.ICP.FRFC.Lon ...
    new_corners_meta.GeoData.ImageCorners.ICP.FRLC.Lon ...
    new_corners_meta.GeoData.ImageCorners.ICP.LRLC.Lon ...
    new_corners_meta.GeoData.ImageCorners.ICP.LRFC.Lon];
for i = 1:4
    [m, corner_ind(i)] = min(sum(bsxfun(@minus, corner_latlons(:,i), new_corner_latlons).^2));
end
if ~isequal(corner_ind, 1:4)
    validation_report = add_val_inc(validation_report, 'Error', ...
        'Corner coordinates ordered incorrectly.', ...
        ['Corner permutation: ' mat2str(corner_ind) ]);
end
% 2.10 GeoData.SCP
ecf1=[SICD_meta.GeoData.SCP.ECF.X ...
    SICD_meta.GeoData.SCP.ECF.Y SICD_meta.GeoData.SCP.ECF.Z];
ecf2=geodetic_to_ecf([SICD_meta.GeoData.SCP.LLH.Lat ...
    SICD_meta.GeoData.SCP.LLH.Lon SICD_meta.GeoData.SCP.LLH.HAE]);
ecf_diff = norm(ecf1(:)-ecf2(:));
if ecf_diff > SCPCOA_TOL
    validation_report = add_val_inc(validation_report, 'Error', ...
        'GeoData.SCP.ECF and GeoData.SCP.LLH not consistent.', ...
        ['GeoData.SCP.ECF - GeoData.SCP.LLH: ' num2str(pos_diff) ' (m)']);
end
% 2.11 ValidData
% Both ImageData.ValidData and GeoData.ValidData are optional, but are
% conditionally required upon each other.  If one exist, the other must
% also.
% 2.11.1
if isfield(SICD_meta.ImageData,'ValidData') && ~isfield(SICD_meta.GeoData,'ValidData')
    validation_report = add_val_inc(validation_report, 'Error', ...
        'ImageData.ValidData/GeoData.ValidData required together.', ...
        'ImageData.ValidData exist but GeoData.ValidData does not.');
end
% 2.11.2
if ~isfield(SICD_meta.ImageData,'ValidData') && isfield(SICD_meta.GeoData,'ValidData')
    validation_report = add_val_inc(validation_report, 'Error', ...
        'ImageData.ValidData/GeoData.ValidData required together.', ...
        'GeoData.ValidData exist but ImageData.ValidData does not.');
end
% As with the corner cordinates, we don't check for precision in the
% agreement between ImageData.ValidData and GeoData.ValidData here, since
% the spec allows from approximate Lat/Lons, and since its not possible to
% know the HAE used to compute Lat/Lon for each corner.  One could
% potentially use the SICD sensor model to check that the image coordinates
% and the resulting ground coordinates from projection to a flat plane (or
% constant HAE) roughly agreed with GeoData.ValidData to within a generous
% tolerance.
if isfield(SICD_meta.ImageData,'ValidData')
    % 2.11.3 In ValidData, first vertex should have (1) minimum row index
    % and (2) minimum column index if 2 vertices exist with equal minimum
    % row index.
    vertex_rows = [SICD_meta.ImageData.ValidData.Vertex.Row];
    vertex_cols = [SICD_meta.ImageData.ValidData.Vertex.Col];
    [min_row, min_vertex_row]=min(vertex_rows);
    if min_vertex_row > 1
        validation_report = add_val_inc(validation_report, 'Error', ...
            'ImageData.ValidData first vertex should have minimum row index.', ...
            ['ImageData.ValidData vertex with minimum row index: ' num2str(min_vertex_row)]);
    else
        [~, min_vertex_col] = min(vertex_cols(vertex_rows==min_row));
        if min_vertex_col > 1
            validation_report = add_val_inc(validation_report, 'Error', ...
                'ImageData.ValidData first vertex should have minimum row/col index.', ...
                ['ImageData.ValidData vertex with minimum row/col index: ' num2str(min_vertex_col)]);
        end
    end
    % 2.11.4 Check that the vertices are ordered clockwise.  ispolycw from
    % the Mapping Toolbox is an easy way to do this if you have the
    % toolbox, but there is code available online to do the same sort of
    % test if we wanted to remove the toolbox dependency.
    % Note: ispolycw assumes Y indexing increases upward, but SICD rows
    % count from top downward, so we must flip sign of vertex_rows.
    if ~isempty(ver('map')) && ~ispolycw(vertex_cols, -vertex_rows)
        validation_report = add_val_inc(validation_report, 'Error', ...
            'ImageData.ValidData vertices must be ordered clockwise.', ...
            'ImageData.ValidData vertices must be ordered clockwise.');
    end
end
    
% 2.12 IFP-specific checks
SCP = [SICD_meta.GeoData.SCP.ECF.X SICD_meta.GeoData.SCP.ECF.Y SICD_meta.GeoData.SCP.ECF.Z];
ruvect = [SICD_meta.Grid.Row.UVectECF.X SICD_meta.Grid.Row.UVectECF.Y SICD_meta.Grid.Row.UVectECF.Z];
cuvect = [SICD_meta.Grid.Col.UVectECF.X SICD_meta.Grid.Col.UVectECF.Y SICD_meta.Grid.Col.UVectECF.Z];
UVECT_TOL = 1e-3;
IFP_POLY_TOL = 1e-5;
switch SICD_meta.ImageFormation.ImageFormAlgo
    case 'RGAZCOMP' % RGAZCOMP case not tested
        % 2.12.1.1
        if ~strcmp(SICD_meta.Grid.ImagePlane,'SLANT')
            validation_report = add_val_inc(validation_report, 'Error', ...
                'RGAZCOMP image formation should result in a SLANT plane image.', ...
                ['Grid.ImagePlane: ' SICD_meta.Grid.ImagePlane ]);
        end
        % 2.12.1.2
        if ~strcmp(SICD_meta.Grid.Type,'RGAZIM')
            validation_report = add_val_inc(validation_report, 'Error', ...
                'RGAZCOMP image formation should result in a RGAZIM grid.', ...
                ['Grid.Type: ' SICD_meta.Grid.Type ]);
        end
        ARP=[SICD_meta.SCPCOA.ARPPos.X SICD_meta.SCPCOA.ARPPos.Y SICD_meta.SCPCOA.ARPPos.Z];
        ARP_v=[SICD_meta.SCPCOA.ARPVel.X SICD_meta.SCPCOA.ARPVel.Y SICD_meta.SCPCOA.ARPVel.Z];
        uRG = (SCP - ARP)/norm(SCP - ARP); % Range unit vector
        left = cross(ARP/norm(ARP),ARP_v/norm(ARP));
        look = sign(left * uRG');
        spn=-look*cross(uRG,ARP_v); spn=spn/norm(spn); % Slant plane unit normal
        % 2.12.1.3
        if ~isfield(SICD_meta,'RgAzComp')
            validation_report = add_val_inc(validation_report, 'Error', ...
                'RGAZCOMP image formation declared, but no RgAzComp metadata given.', ...
                'No RgAzComp field in this SICD.');
        else
            % 2.12.1.4
            derived_AzSF = -look * sind(SICD_meta.SCPCOA.DopplerConeAng)/...
                SICD_meta.SCPCOA.SlantRange;
            if abs(SICD_meta.RgAzComp.AzSF - derived_AzSF) > 1e-6
                validation_report = add_val_inc(validation_report, 'Error', ...
                    'RgAzComp fields inconsistent.', ...
                    ['RgAzComp.AzSF: ' num2str(SICD_meta.RgAzComp.AzSF) ...
                    ' , Derived RgAzComp.AzSF: ' num2str(derived_AzSF) ]);
            end
            % 2.12.1.5
            if isfield(SICD_meta.Timeline,'IPP') && ...
                    numel(SICD_meta.Timeline.IPP.Set)==1
                krg_coa = SICD_meta.Grid.Row.KCtr;
                if isfield(SICD_meta.Grid.Row,'DeltaKCOAPoly')
                    krg_coa = krg_coa + SICD_meta.Grid.Row.DeltaKCOAPoly;
                end
                st_rate_coa = polyval(polyder(SICD_meta.Timeline.IPP.Set.IPPPoly(end:-1:1)),SICD_meta.SCPCOA.SCPTime);
                delta_kaz_per_delta_v = look * krg_coa * ... % Scale factor described in SICD spec
                    (norm(ARP_v) * sind(SICD_meta.SCPCOA.DopplerConeAng) / SICD_meta.SCPCOA.SlantRange) / ...
                    st_rate_coa;
                derived_KazPoly = delta_kaz_per_delta_v * ...
                    SICD_meta.Timeline.IPP.Set.IPPPoly;
                if norm(SICD_meta.RgAzComp.KazPoly - derived_KazPoly) > 1e-3
                    validation_report = add_val_inc(validation_report, 'Error', ...
                        'RgAzComp fields inconsistent.', ...
                        ['RgAzComp.KazPoly: ' mat2str(SICD_meta.RgAzComp.KazPoly) ...
                        ' , Derived RgAzComp.KazPoly: ' mat2str(derived_KazPoly) ]);
                end
            end
        end
        % 2.12.1.6
        drvect = uRG;
        if norm(drvect(:) - ruvect(:)) > UVECT_TOL
            validation_report = add_val_inc(validation_report, 'Error', ...
                'UVect fields inconsistent.', ...
                ['Grid.Row.UVectECF: ' mat2str(ruvect) ...
                ', Derived Grid.Row.UVectECF: ' mat2str(drvect)]);
        end
        % 2.12.1.7
        dcvect = cross(spn,uRG);
        if norm(dcvect(:) - cuvect(:)) > UVECT_TOL
            validation_report = add_val_inc(validation_report, 'Error', ...
                'UVect fields inconsistent.', ...
                ['Grid.Col.UVectECF: ' mat2str(cuvect) ...
                ', Derived Grid.Col.UVectECF: ' mat2str(dcvect)]);
        end
        % 2.12.1.8
        if isfield(SICD_meta.Grid.Col,'DeltaKCOAPoly') && ...
                isscalar(SICD_meta.Grid.Col.DeltaKCOAPoly) && ...
                (abs(SICD_meta.Grid.Col.KCtr - (-SICD_meta.Grid.Col.DeltaKCOAPoly)) > ...
                eps(SICD_meta.Grid.Col.KCtr))
            validation_report = add_val_inc(validation_report, 'Error', ...
                'Grid.Col.KCtr must be equal to -Grid.Col.DeltaKCOAPoly for RGAZCOMP data.', ...
                ['Grid.Col.KCtr: ' num2str(SICD_meta.Grid.Col.KCtr) ...
                ', Grid.Col.DeltaKCOAPoly: ' num2str(SICD_meta.Grid.Col.DeltaKCOAPoly)]);
        elseif SICD_meta.Grid.Col.KCtr ~= 0
            validation_report = add_val_inc(validation_report, 'Error', ...
                'Grid.Col.KCtr must be zero for RGAZCOMP data with no Grid.Col.DeltaKCOAPoly field.', ...
                ['Grid.Col.KCtr = ' num2str(SICD_meta.Grid.Col.KCtr)]);
        end
        % 2.12.1.9
        if ~isfield(SICD_meta.RadarCollection,'RefFreqIndex') 
            fc_proc = (SICD_meta.ImageFormation.TxFrequencyProc.MinProc + ...
                SICD_meta.ImageFormation.TxFrequencyProc.MaxProc)/2;
            k_f_c = fc_proc * (2/SPEED_OF_LIGHT);
            if isfield(SICD_meta.Grid.Row,'DeltaKCOAPoly') && ...
                    isscalar(SICD_meta.Grid.Col.DeltaKCOAPoly) && ...
                    abs(SICD_meta.Grid.Row.KCtr - (k_f_c - SICD_meta.Grid.Row.DeltaKCOAPoly)) > eps(k_f_c)
                validation_report = add_val_inc(validation_report, 'Error', ...
                    WF_INCONSISTENT_STR, ...
                    ['Grid.Row.KCtr: ' num2str(SICD_meta.Grid.Row.KCtr) ...
                    ', Center frequency * 2/c - Grid.Row.DeltaKCOAPoly: ' ...
                    num2str(k_f_c - SICD_meta.Grid.Row.DeltaKCOAPoly)]);
            elseif abs(SICD_meta.Grid.Row.KCtr - k_f_c) > eps(k_f_c)
                validation_report = add_val_inc(validation_report, 'Error', ...
                    WF_INCONSISTENT_STR, ...
                    ['Grid.Row.KCtr: ' num2str(SICD_meta.Grid.Col.KCtr) ...
                    ', Center frequency * 2/c: ' num2str(k_f_c)]);
            end
        end
        % 2.12.1.10
        if isfield(SICD_meta.Grid.Col,'DeltaKCOAPoly') && ...
                any(SICD_meta.Grid.Col.DeltaKCOAPoly(2:end))
            validation_report = add_val_inc(validation_report, 'Error', ...
                'Grid.Col.DeltaKCOAPoly must be a single value for RGAZCOMP data.', ...
                ['Grid.Col.DeltaKCOAPoly = ' mat2str(SICD_meta.Grid.Col.DeltaKCOAPoly)]);
        end
        % 2.12.1.11
        if isfield(SICD_meta.Grid.Row,'DeltaKCOAPoly') && ...
                any(SICD_meta.Grid.Row.DeltaKCOAPoly(2:end))
            validation_report = add_val_inc(validation_report, 'Error', ...
                'Grid.Row.DeltaKCOAPoly must be a single value for RGAZCOMP data.', ...
                ['Grid.Row.DeltaKCOAPoly = ' mat2str(SICD_meta.Grid.Row.DeltaKCOAPoly)]);
        end
    case 'PFA'
        % 2.12.2.1
        if ~strcmp(SICD_meta.Grid.Type,'RGAZIM')
            validation_report = add_val_inc(validation_report, 'Error', ...
                'PFA image formation should result in a RGAZIM grid.', ...
                ['Grid.Type: ' SICD_meta.Grid.Type ]);
        end
        % 2.12.2.2
        if ~isfield(SICD_meta,'PFA')
            validation_report = add_val_inc(validation_report, 'Error', ...
                'PFA image formation declared, but no PFA metadata given.', ...
                'No PFA field in this SICD.');
        else
            % 2.12.2.3
            if (SICD_meta.PFA.PolarAngRefTime - SICD_meta.SCPCOA.SCPTime) > eps
                validation_report = add_val_inc(validation_report, 'Warning', ...
                    'Polar angle reference time and center of aperture time for scene center are usually the same.', ...
                    ['PFA.PolarAngRefTime: ' ...
                    num2str(SICD_meta.PFA.PolarAngRefTime) ...
                    ' , SCPCOA.SCPTime: ' ...
                    num2str(SICD_meta.SCPCOA.SCPTime) ]);
            end
            % Check PFA polynomials
            num_pfa_coefs = numel(SICD_meta.PFA.PolarAngPoly);
            times = linspace(0, SICD_meta.Timeline.CollectDuration, num_pfa_coefs+1).';
            pos = [polyval(SICD_meta.Position.ARPPoly.X(end:-1:1), times) ...
                polyval(SICD_meta.Position.ARPPoly.Y(end:-1:1), times) ...
                polyval(SICD_meta.Position.ARPPoly.Z(end:-1:1), times)];
            pol_ref_pos = [polyval(SICD_meta.Position.ARPPoly.X(end:-1:1), SICD_meta.PFA.PolarAngRefTime) ...
                polyval(SICD_meta.Position.ARPPoly.Y(end:-1:1), SICD_meta.PFA.PolarAngRefTime) ...
                polyval(SICD_meta.Position.ARPPoly.Z(end:-1:1), SICD_meta.PFA.PolarAngRefTime)];
            ipn = [SICD_meta.PFA.IPN.X SICD_meta.PFA.IPN.Y SICD_meta.PFA.IPN.Z];
            fpn = [SICD_meta.PFA.FPN.X SICD_meta.PFA.FPN.Y SICD_meta.PFA.FPN.Z];
            [k_a, k_sf] = pfa_polar_coords(pos, SCP, pol_ref_pos, ipn, fpn);
            old_state = warning('off','MATLAB:polyfit:RepeatedPointsOrRescale');
            PolarAngPoly = fliplr( polyfit(times, k_a, num_pfa_coefs - 1) ).';
            SpatialFreqSFPoly = fliplr( polyfit(k_a, k_sf, num_pfa_coefs - 1) ).';
            warning(old_state);
            % 2.12.2.4
            if max(abs(polyval(SICD_meta.PFA.PolarAngPoly(end:-1:1),times)-k_a))>IFP_POLY_TOL
                validation_report = add_val_inc(validation_report, 'Error', ...
                    'PFA fields inconsistent.', ...
                    ['PFA.PolarAngPoly: ' mat2str(SICD_meta.PFA.PolarAngPoly) ...
                    ' , Derived PolarAngPoly: ' mat2str(PolarAngPoly)]);
            end
            % 2.12.2.5
            if max(abs(polyval(SICD_meta.PFA.SpatialFreqSFPoly(end:-1:1),k_a)-k_sf))>IFP_POLY_TOL
                validation_report = add_val_inc(validation_report, 'Error', ...
                    'PFA fields inconsistent.', ...
                    ['PFA.SpatialFreqSFPoly: ' mat2str(SICD_meta.PFA.SpatialFreqSFPoly) ...
                    ' , Derived SpatialFreqSFPoly: ' mat2str(SpatialFreqSFPoly)]);
            end
            % 2.12.2.6
            % Row.UVect should be the range vector at zero polar angle time
            % projected into IPN
            % Projection of a point along a given direction to a plane is
            % just the intersection of the line defined by that point (l0)
            % and direction (l) and the plane defined by a point in the
            % plane (p0) and the normal (p):
            % l0 + ((l0 - p0).p/(l.p))*l
            % where . represents the dot product.
            d = (bsxfun(@minus, SCP, pol_ref_pos) * ipn(:)) ./...
                (fpn * ipn(:)); % Distance from point to plane in line_direction
            ref_pos_ipn = pol_ref_pos + (d * fpn);
            srv=(ref_pos_ipn-SCP).'; % slant range vector
            srv_unit=srv/norm(srv); % slant range unit vector
            drvect = -srv_unit;
            if norm(drvect(:) - ruvect(:)) > UVECT_TOL
                validation_report = add_val_inc(validation_report, 'Error', ...
                    'UVect fields inconsistent.', ...
                    ['Grid.Row.UVectECF: ' mat2str(ruvect) ...
                    ', Derived Grid.Row.UVectECF: ' mat2str(drvect)]);
            end
            % 2.12.2.7
            dcvect = cross(srv_unit,ipn);
            if norm(dcvect(:) - cuvect(:)) > UVECT_TOL
                validation_report = add_val_inc(validation_report, 'Error', ...
                    'UVect fields inconsistent.', ...
                    ['Grid.Col.UVectECF: ' mat2str(cuvect) ...
                    ', Derived Grid.Col.UVectECF: ' mat2str(dcvect)]);
            end
            % 2.12.2.8
            if isfield(SICD_meta.PFA,'STDeskew') && ...
                    SICD_meta.PFA.STDeskew.Applied && ...
                    ~any(SICD_meta.Grid.TimeCOAPoly(2:end))
                validation_report = add_val_inc(validation_report, 'Error', ...
                    'STDeskew only appropriate for data with varying Center of Aperture times.', ...
                    ['TimeCOAPoly: ' mat2str(SICD_meta.Grid.TimeCOAPoly) ...
                    ', PFA.STDeskew.Appled: True']);
            end
            % 2.12.2.9
            % Row.DeltaKCOAPoly should be derivative of STDeskew (if it exist)
            if isfield(SICD_meta.PFA,'STDeskew') && ...
                    SICD_meta.PFA.STDeskew.Applied && ...
                    isfield(SICD_meta.Grid.Row, 'DeltaKCOAPoly')
                STDSPhasePoly = SICD_meta.PFA.STDeskew.STDSPhasePoly;
                DeltaKCOAPoly = SICD_meta.Grid.Row.DeltaKCOAPoly;
                new_DeltaKCOAPoly = STDSPhasePoly(2:end,:) .* ...
                    repmat((1:(size(STDSPhasePoly,1)-1)).',1,size(STDSPhasePoly,2));
                % Sometime Grid.Row.DeltaKCOAPoly and new_DeltaKCOAPoly
                % could vary in size.  Difference should be zeropad only.
                bigger_size = max(size(new_DeltaKCOAPoly), size(DeltaKCOAPoly));
                padded_DeltaKCOAPoly = zeros(bigger_size);
                padded_derived_DeltaKCOAPoly = zeros(bigger_size);
                padded_DeltaKCOAPoly(1:size(DeltaKCOAPoly,1), 1:size(DeltaKCOAPoly,2)) = DeltaKCOAPoly;
                padded_derived_DeltaKCOAPoly(1:size(new_DeltaKCOAPoly,1), 1:size(new_DeltaKCOAPoly,2)) = new_DeltaKCOAPoly;
                if max(abs(padded_DeltaKCOAPoly(:) - padded_derived_DeltaKCOAPoly(:))) > IFP_POLY_TOL
                    validation_report = add_val_inc(validation_report, 'Error', ...
                        'Slow-time deskew polynomial inconsistent.', ...
                        ['Grid.Row.DeltaKCOAPoly: ' mat2str(SICD_meta.Grid.Row.DeltaKCOAPoly) ...
                        ', Derived Grid.Row.DeltaKCOAPoly: ' mat2str(new_DeltaKCOAPoly)]);
                end
            end
            % 2.12.2.10
            % Make sure Row.KCtr is consistent with processed RF frequency bandwidth
            if ~isfield(SICD_meta.RadarCollection,'RefFreqIndex')
                fc_proc = ((SICD_meta.ImageFormation.TxFrequencyProc.MinProc + ...
                    SICD_meta.ImageFormation.TxFrequencyProc.MaxProc)/2);
                % PFA.SpatialFreqSFPoly affects Row.KCtr.
                kap_ctr = fc_proc*SICD_meta.PFA.SpatialFreqSFPoly(1)*2/SPEED_OF_LIGHT;
                % PFA inscription could cause kap_ctr and Row.KCtr to be somewhat different
                theta = atan((SICD_meta.Grid.Col.ImpRespBW/2) / SICD_meta.Grid.Row.KCtr); % Aperture angle
                kctr_tol = 1 - cos(theta); % Difference between Krg and Kap (Krg = cos(theta)*Kap)
                kctr_tol = max(0.01, kctr_tol); % .01 should be plenty of tolerance for precision issues at small angles
                if abs((SICD_meta.Grid.Row.KCtr/kap_ctr) - 1) > kctr_tol
                    validation_report = add_val_inc(validation_report, 'Error', ...
                        WF_INCONSISTENT_STR, ...
                        ['Grid.Row.KCtr: ' num2str(SICD_meta.Grid.Row.KCtr) ...
                        ', Derived Kap_Ctr: ' num2str(kap_ctr)]);
                end
            end
            % 2.12.2.11 Check origin of PolarAngPoly for consistency
            polar_ang_ref = polyval(SICD_meta.PFA.PolarAngPoly(end:-1:1), ...\
                SICD_meta.PFA.PolarAngRefTime); % Should be zero
            if abs(polar_ang_ref)>IFP_POLY_TOL
                validation_report = add_val_inc(validation_report, 'Error', ...
                    'PFA fields inconsistent.', ...
                    ['PFA.PolarAngPoly evaluated at PFA.PolarAngRefTime: ' ...
                    num2str(polar_ang_ref)]);
            end
            % 2.12.2.12 Check that image plane normal points away from center of earth
            if dot(ipn,SCP) < 0
                validation_report = add_val_inc(validation_report, 'Error', ...
                    'Image formation plane unit normal must point away from center of earth.', ...
                    ['PFA.IPN: ' num2str(ipn) ]);
            end
            % 2.12.2.13 Check whether image plane is roughly close to claimed Grid.ImagePlane
            if strcmpi(SICD_meta.Grid.ImagePlane,'SLANT')
                % Cut-and-paste from RGAZCOMP section
                ARP=[SICD_meta.SCPCOA.ARPPos.X SICD_meta.SCPCOA.ARPPos.Y SICD_meta.SCPCOA.ARPPos.Z];
                ARP_v=[SICD_meta.SCPCOA.ARPVel.X SICD_meta.SCPCOA.ARPVel.Y SICD_meta.SCPCOA.ARPVel.Z];
                uRG = (SCP - ARP)/norm(SCP - ARP); % Range unit vector
                left = cross(ARP/norm(ARP),ARP_v/norm(ARP));
                look = sign(left * uRG');
                spn=-look*cross(uRG,ARP_v); spn=spn/norm(spn); % Slant plane unit normal
                if acosd(dot(ipn/norm(ipn),spn)) > 2
                    validation_report = add_val_inc(validation_report, 'Error', ...
                        'Image plane claimed to be ''SLANT'', but not close to instantaneous slant plane at COA.', ...
                        ['PFA.IPN: ' num2str(ipn) ...
                        ', COA Slant: ' num2str(spn)]);
                end
            elseif strcmpi(SICD_meta.Grid.ImagePlane,'GROUND')
                if acosd(dot(ipn/norm(ipn),wgs_84_norm(SCP))) > 5
                    validation_report = add_val_inc(validation_report, 'Error', ...
                        'Image plane claimed to be ''GROUND'', but not close to tangent to ellipsoid.', ...
                        ['PFA.IPN: ' num2str(ipn) ...
                        ', Normal to WGS_84: ' num2str(reshape(wgs_84_norm(SCP),1,[]))]);
                end
            end
            % 2.12.2.14 Check that focus plane normal points away from center of earth
            if dot(fpn,SCP) < 0
                validation_report = add_val_inc(validation_report, 'Error', ...
                    'Focus plane unit normal must point away from center of earth.', ...
                    ['PFA.FPN: ' num2str(fpn) ]);
            end
            % 2.12.2.15 Check whether focus plane is close to (within 5
            % degrees) tangent to elliposid
            if acosd(dot(fpn/norm(fpn),wgs_84_norm(SCP))) > 5
                validation_report = add_val_inc(validation_report, 'Warning', ...
                    'Focus plane unit normal generally tangent to ellipsoid.', ...
                    ['PFA.FPN: ' num2str(fpn) ...
                    ', Normal to WGS_84: ' num2str(reshape(wgs_84_norm(SCP),1,[]))]);
            end
            % 2.12.2.16 Check for polar angle consistency
            if strcmp(SICD_meta.CollectionInfo.RadarMode.ModeType,'SPOTLIGHT')
                % In spotlight mode, frequency support covers total polar angle span
                % From polar angle polynomial (previously checked against ARPPoly)
                polar_angle_bounds1 = sort(polyval(SICD_meta.PFA.PolarAngPoly(end:-1:1), ...
                    [SICD_meta.ImageFormation.TStartProc, SICD_meta.ImageFormation.TEndProc]));
                % From frequency support (previously checked against KCtr/ImpRespBW)
                polar_angle_bounds2 = atan([SICD_meta.PFA.Kaz1 SICD_meta.PFA.Kaz2] / ...
                    SICD_meta.PFA.Krg1);
                if any(abs(rad2deg(polar_angle_bounds1-polar_angle_bounds2))>.01)
                    % This equality will only be exactly true for an
                    % inscribed rectangular aperture.
                    validation_report = add_val_inc(validation_report, 'Error', ...
                        'Polar angle bounds inconsistent.', ...
                        ['From PFA.PolarAngPoly (rad): ' num2str(polar_angle_bounds1) ...
                        ', From PFA.Kaz/Krg (rad): ' num2str(polar_angle_bounds2)]);
                end
            end
        end
    case 'RMA'
        % 2.12.3.1
        if ~isfield(SICD_meta,'RMA')
            validation_report = add_val_inc(validation_report, 'Error', ...
                'RMA image formation declared, but no RMA metadata given.', ...
                'No RMA field in this SICD.');
        else
            switch SICD_meta.RMA.ImageType
                case 'RMAT'
                    % 2.12.3.2.1
                    if ~strcmp(SICD_meta.Grid.Type,'XCTYAT')
                        validation_report = add_val_inc(validation_report, 'Error', ...
                            'RMA/RMAT image formation should result in a XCTYAT grid.', ...
                            ['RMA.ImageType: RMAT, Grid.Type: ' SICD_meta.Grid.Type ]);
                    end
                    % 2.12.3.2.2
                    if ~isfield(SICD_meta.RMA,'RMAT')
                        validation_report = add_val_inc(validation_report, 'Error', ...
                            'RMA/RMAT image formation declared, but no RMAT metadata given.', ...
                            'No RMAT field in this SICD.');
                    else
                        PosRef=[SICD_meta.RMA.RMAT.PosRef.X SICD_meta.RMA.RMAT.PosRef.Y SICD_meta.RMA.RMAT.PosRef.Z];
                        VelRef=[SICD_meta.RMA.RMAT.VelRef.X SICD_meta.RMA.RMAT.VelRef.Y SICD_meta.RMA.RMAT.VelRef.Z];
                        uLOS = (SCP - PosRef)/norm(SCP - PosRef); % Range unit vector
                        left = cross(PosRef/norm(PosRef),VelRef/norm(VelRef));
                        look = sign(left * uLOS');
                        uYAT = -look*VelRef/norm(VelRef); % Along track unit vector
                        spn = cross(uLOS,uYAT); spn=spn/norm(spn); % Reference slant plane normal
                        uXCT = cross(uYAT, spn); % Cross track unit vector
                        dcaRef = acosd((VelRef/norm(VelRef))*(uLOS.'));
                        % 2.12.3.2.3
                        if norm(ruvect(:) - uXCT(:)) > UVECT_TOL
                            validation_report = add_val_inc(validation_report, 'Error', ...
                                'UVect fields inconsistent.', ...
                                ['Grid.Row.UVectECF: ' mat2str(ruvect) ...
                                ', Derived Grid.Row.UVectECF: ' mat2str(uXCT)]);
                        end
                        % 2.12.3.2.4
                        if norm(cuvect(:) - uYAT(:)) > UVECT_TOL
                            validation_report = add_val_inc(validation_report, 'Error', ...
                                'UVect fields inconsistent.', ...
                                ['Grid.Col.UVectECF: ' mat2str(cuvect) ...
                                ', Derived Grid.Col.UVectECF: ' mat2str(uYAT)]);
                        end
                        % 2.12.3.2.5
                        if abs(dcaRef - SICD_meta.RMA.RMAT.DopConeAngRef) > 1e-6
                            validation_report = add_val_inc(validation_report, 'Error', ...
                                'RMA fields inconsistent.', ...
                                ['RMA.RMAT.DopConeAngRef: ' num2str(SICD_meta.RMA.RMAT.DopConeAngRef) ...
                                ', Derived RMA.RMAT.DopConeAngRef: ' num2str(dcaRef)]);
                        end
                    end
                    if ~isfield(SICD_meta.RadarCollection,'RefFreqIndex') 
                        fc_proc = (SICD_meta.ImageFormation.TxFrequencyProc.MinProc + ...
                            SICD_meta.ImageFormation.TxFrequencyProc.MaxProc)/2;
                        k_f_c = fc_proc * (2/SPEED_OF_LIGHT);
                        % 2.12.3.2.6
                        if abs((k_f_c * sind(SICD_meta.RMA.RMAT.DopConeAngRef) / ...
                                SICD_meta.Grid.Row.KCtr) - 1) > WF_TOL
                            validation_report = add_val_inc(validation_report, 'Warning', ...
                                WF_INCONSISTENT_STR, ...
                                ['Grid.Row.KCtr = ' num2str(SICD_meta.Grid.Row.KCtr) ...
                                ', Derived Grid.Row.KCtr: ' num2str(...
                                k_f_c * sind(SICD_meta.RMA.RMAT.DopConeAngRef))]);
                        end
                        % 2.12.3.2.7
                        if abs((k_f_c * cosd(SICD_meta.RMA.RMAT.DopConeAngRef) / ...
                                SICD_meta.Grid.Col.KCtr) - 1) > WF_TOL
                            validation_report = add_val_inc(validation_report, 'Warning', ...
                                WF_INCONSISTENT_STR, ...
                                ['Grid.Col.KCtr = ' num2str(SICD_meta.Grid.Col.KCtr) ...
                                ', Derived Grid.Col.KCtr: ' num2str(...
                                k_f_c * cosd(SICD_meta.RMA.RMAT.DopConeAngRef))]);
                        end
                    end
                case 'RMCR'
                    % 2.12.3.3.1
                    if ~strcmp(SICD_meta.Grid.Type,'XRGYCR')
                        validation_report = add_val_inc(validation_report, 'Error', ...
                            'RMA/RMCR image formation should result in a XRGYCR grid.', ...
                            ['RMA.ImageType: RMCR, Grid.Type: ' SICD_meta.Grid.Type ]);
                    end
                    % 2.12.3.3.2
                    if ~isfield(SICD_meta.RMA,'RMCR')
                        validation_report = add_val_inc(validation_report, 'Error', ...
                            'RMA/RMCR image formation declared, but no RMCR metadata given.', ...
                            'No RMCR field in this SICD.');
                    else
                        PosRef = [SICD_meta.RMA.RMCR.PosRef.X ...
                            SICD_meta.RMA.RMCR.PosRef.Y ...
                            SICD_meta.RMA.RMCR.PosRef.Z];
                        VelRef = [SICD_meta.RMA.RMCR.VelRef.X ...
                            SICD_meta.RMA.RMCR.VelRef.Y ...
                             SICD_meta.RMA.RMCR.VelRef.Z];
                        uXRG = (SCP - PosRef)/norm(SCP - PosRef); % Range unit vector
                        left = cross(PosRef/norm(PosRef),VelRef/norm(VelRef));
                        look = sign(left * uXRG');
                        spn = look * cross(VelRef / norm(VelRef), uXRG);
                        spn=spn/norm(spn); % Reference slant plane normal
                        uYCR = cross(spn, uXRG); % Cross range unit vector
                        dcaRef = acosd((VelRef/norm(VelRef))*(uXRG.'));
                         % 2.12.3.3.3
                        if norm(ruvect(:) - uXRG(:)) > UVECT_TOL
                            validation_report = add_val_inc(validation_report, 'Error', ...
                                'UVect fields inconsistent.', ...
                                ['Grid.Row.UVectECF: ' mat2str(ruvect) ...
                                ', Derived Grid.Row.UVectECF: ' mat2str(uXRG)]);
                        end
                         % 2.12.3.3.4
                        if norm(cuvect(:) - uYCR(:)) > UVECT_TOL
                            validation_report = add_val_inc(validation_report, 'Error', ...
                                'UVect fields inconsistent.', ...
                                ['Grid.Col.UVectECF: ' mat2str(cuvect) ...
                                ', Derived Grid.Col.UVectECF: ' mat2str(uYCR)]);
                        end
                         % 2.12.3.3.5
                        if abs(dcaRef - SICD_meta.RMA.RMCR.DopConeAngRef) > 1e-6
                            validation_report = add_val_inc(validation_report, 'Error', ...
                                'RMA fields inconsistent.', ...
                                ['RMA.RMCR.DopConeAngRef: ' num2str(SICD_meta.RMA.RMCR.DopConeAngRef) ...
                                ', Derived RMA.RMCR.DopConeAngRef: ' num2str(dcaRef)]);
                        end
                    end
                    % 2.12.3.3.6
                    if SICD_meta.Grid.Col.KCtr ~= 0
                        validation_report = add_val_inc(validation_report, 'Error', ...
                            'Grid.Col.KCtr must be zero for RMA/RMCR data.', ...
                            ['Grid.Col.KCtr = ' num2str(SICD_meta.Grid.Col.KCtr)]);
                    end
                    % 2.12.3.3.7
                    if ~isfield(SICD_meta.RadarCollection,'RefFreqIndex') 
                        fc_proc = (SICD_meta.ImageFormation.TxFrequencyProc.MinProc + ...
                            SICD_meta.ImageFormation.TxFrequencyProc.MaxProc)/2;
                        k_f_c = fc_proc * (2/SPEED_OF_LIGHT);
                        if abs((SICD_meta.Grid.Row.KCtr / k_f_c) - 1) > WF_TOL
                            validation_report = add_val_inc(validation_report, 'Warning', ...
                                WF_INCONSISTENT_STR, ...
                                ['Grid.Row.KCtr = ' num2str(SICD_meta.Grid.Row.KCtr) ...
                                ', Center frequency * 2/c: ' num2str(k_f_c)]);
                        end
                    end
                case 'INCA'
                    % 2.12.3.4.1
                    if ~strcmp(SICD_meta.Grid.Type,'RGZERO')
                        validation_report = add_val_inc(validation_report, 'Error', ...
                            'RMA/INCA image formation should result in a RGZERO grid.', ...
                            ['Grid.Type: ' SICD_meta.Grid.Type ]);
                    end
                    % 2.12.3.4.2
                    if ~isfield(SICD_meta.RMA,'INCA')
                        validation_report = add_val_inc(validation_report, 'Error', ...
                            'RMA/INCA image formation declared, but no INCA metadata given.', ...
                            'No INCA field in this SICD.');
                    else
                        % 2.12.3.4.3
                        if strcmp(SICD_meta.CollectionInfo.RadarMode.ModeType,'SPOTLIGHT') && ...
                            (isfield(SICD_meta.RMA.INCA, 'DopCentroidPoly') || ...
                            isfield(SICD_meta.RMA.INCA, 'DopCentroidCOA'))
                            validation_report = add_val_inc(validation_report, 'Error', ...
                                'RMA.INCA fields inconsistent.', ...
                                'RMA.INCA.DopCentroidPoly/DopCentroidCOA not used for SPOTLIGHT collection.');
                        % 2.12.3.4.4
                        elseif ~strcmp(SICD_meta.CollectionInfo.RadarMode.ModeType,'SPOTLIGHT')
                            if (~isfield(SICD_meta.RMA.INCA, 'DopCentroidPoly') || ...
                                    ~isfield(SICD_meta.RMA.INCA, 'DopCentroidCOA'))
                                validation_report = add_val_inc(validation_report, 'Error', ...
                                    'RMA.INCA fields inconsistent.', ...
                                    'RMA.INCA.DopCentroidPoly/DopCentroidCOA required for non-SPOTLIGHT collection.');
                            % 2.12.3.4.5
                            elseif isfield(SICD_meta.Grid.Col,'DeltaKCOAPoly') && ...
                                    SICD_meta.RMA.INCA.DopCentroidCOA && ...
                                    (norm(SICD_meta.Grid.Col.DeltaKCOAPoly - ...
                                    (SICD_meta.RMA.INCA.DopCentroidPoly * ...
                                    SICD_meta.RMA.INCA.TimeCAPoly(2))) > IFP_POLY_TOL)
                                validation_report = add_val_inc(validation_report, 'Error', ...
                                    'RMA.INCA fields inconsistent.', ...
                                    ['Grid.Col.DeltaKCOAPoly: ' mat2str(SICD_meta.Grid.Col.DeltaKCOAPoly) ...
                                    ', RMA.INCA.DopCentroidPoly*RMA.INCA.TimeCAPoly(2): ' ...
                                    mat2str(SICD_meta.RMA.INCA.DopCentroidPoly*...
                                    SICD_meta.RMA.INCA.TimeCAPoly(2))]);
                            end
                        end
                    end
                    % INCA UVects are defined from closest approach
                    % position/velocity, not center of aperture
                    ca_pos = [polyval(SICD_meta.Position.ARPPoly.X(end:-1:1), ...
                        SICD_meta.RMA.INCA.TimeCAPoly(1)) ...
                        polyval(SICD_meta.Position.ARPPoly.Y(end:-1:1), ...
                        SICD_meta.RMA.INCA.TimeCAPoly(1)) ...
                        polyval(SICD_meta.Position.ARPPoly.Z(end:-1:1), ...
                        SICD_meta.RMA.INCA.TimeCAPoly(1))];
                    ca_vel = [polyval(polyder(SICD_meta.Position.ARPPoly.X(end:-1:1)), ...
                        SICD_meta.RMA.INCA.TimeCAPoly(1)) ...
                        polyval(polyder(SICD_meta.Position.ARPPoly.Y(end:-1:1)), ...
                        SICD_meta.RMA.INCA.TimeCAPoly(1)) ...
                        polyval(polyder(SICD_meta.Position.ARPPoly.Z(end:-1:1)), ...
                        SICD_meta.RMA.INCA.TimeCAPoly(1))];
                    uRG = (SCP - ca_pos)/norm(SCP - ca_pos); % Range unit vector
                    left = cross(ca_pos/norm(ca_pos),ca_vel/norm(ca_pos));
                    look = sign(left * uRG');
                    spn=-look*cross(uRG,ca_vel); spn=spn/norm(spn); % Slant plane unit normal
                    uAZ = cross(spn,uRG);
                    % 2.12.3.4.6
                    if norm(uRG(:) - ruvect(:)) > UVECT_TOL
                        validation_report = add_val_inc(validation_report, 'Error', ...
                            'UVect fields inconsistent.', ...
                            ['Grid.Row.UVectECF: ' mat2str(ruvect) ...
                            ', Derived Grid.Row.UVectECF: ' mat2str(uRG)]);
                    end
                    % 2.12.3.4.7
                    if norm(uAZ(:) - cuvect(:)) > UVECT_TOL
                        validation_report = add_val_inc(validation_report, 'Error', ...
                            'UVect fields inconsistent.', ...
                            ['Grid.Col.UVectECF: ' mat2str(cuvect) ...
                            ', Derived Grid.Col.UVectECF: ' mat2str(uAZ)]);
                    end
                    % 2.12.3.4.8
                    if SICD_meta.Grid.Col.KCtr ~= 0
                        validation_report = add_val_inc(validation_report, 'Error', ...
                            'Grid.Col.KCtr must be zero for RMA/INCA data.', ...
                            ['Grid.Col.KCtr = ' num2str(SICD_meta.Grid.Col.KCtr)]);
                    end
                    % 2.12.3.4.9
                    fc = ((SICD_meta.RadarCollection.TxFrequency.Min + ...
                            SICD_meta.RadarCollection.TxFrequency.Max)/2);
                    if abs((SICD_meta.RMA.INCA.FreqZero / fc) - 1) > WF_TOL
                        validation_report = add_val_inc(validation_report, 'Warning', ...
                            'RMA.INCA.FreqZero is typically the center transmit frequency.', ...
                            ['RMA.INCA.FreqZero = ' num2str(SICD_meta.RMA.INCA.FreqZero) ...
                            ', Center transmit frequency: ' num2str(fc)]);
                    end
                    % 2.12.3.4.10
                    if abs(norm(ca_pos-SCP) - SICD_meta.RMA.INCA.R_CA_SCP) > SCPCOA_TOL
                        validation_report = add_val_inc(validation_report, 'Error', ...
                            'RMA.INCA fields inconsistent.', ...
                            ['RMA.INCA.R_CA_SCP = ' num2str(SICD_meta.RMA.INCA.R_CA_SCP) ...
                            ', Derived RMA.INCA.R_CA_SCP: ' num2str(norm(ca_pos-SCP))]);
                        
                    end
                    % 2.12.3.4.11
                    if ~isfield(SICD_meta.RadarCollection,'RefFreqIndex') && ...
                            abs((SICD_meta.Grid.Row.KCtr) - ...
                            SICD_meta.RMA.INCA.FreqZero*2/SPEED_OF_LIGHT) > ...
                            eps(SICD_meta.Grid.Row.KCtr);
                        validation_report = add_val_inc(validation_report, 'Error', ...
                            WF_INCONSISTENT_STR, ...
                            ['RMA.INCA.FreqZero * 2/c: ' num2str(SICD_meta.RMA.INCA.FreqZero*2/SPEED_OF_LIGHT), ...
                            ', Grid.Row.KCtr: ' num2str(SICD_meta.Grid.Row.KCtr)]);
                    end
                    % 2.12.3.4.12
                    if abs((1/(norm(ca_vel)*abs(SICD_meta.RMA.INCA.TimeCAPoly(2)))) - ...
                            SICD_meta.RMA.INCA.DRateSFPoly(1,1)) > 0.001
                        validation_report = add_val_inc(validation_report, 'Error', ...
                            'RMA.INCA fields inconsistent.', ...
                            ['RMA.INCA.DRateSFPoly(1,1) = ' ...
                            num2str(SICD_meta.RMA.INCA.DRateSFPoly(1,1)) ...
                            ', Derived RMA.INCA.DRateSFPoly(1,1): ' ...
                            num2str(1/(norm(ca_vel)*SICD_meta.RMA.INCA.TimeCAPoly(2)))]);
                        
                    end                    
            end
        end
    case 'OTHER'
        % 2.12.3
        validation_report = add_val_inc(validation_report, 'Warning', ...
            'Image formation not fully defined.', ...
            'SICD_meta.ImageFormation.ImageFormAlgo = ''OTHER''');
end
% 2.12.4
if ~(SICD_meta.ImageFormation.TStartProc<SICD_meta.ImageFormation.TEndProc)
    validation_report = add_val_inc(validation_report, 'Error', ...
        'Image formation slow time bounds not valid.', ...
        ['SICD_meta.ImageFormation.TStartProc = ' ...
        num2str(SICD_meta.ImageFormation.TStartProc) ...
        ', SICD_meta.ImageFormation.TEndProc: ' ...
        num2str(SICD_meta.ImageFormation.TEndProc)]);
end
% 2.12.5
if ~ismember(SICD_meta.ImageFormation.TxRcvPolarizationProc, ...
        {SICD_meta.RadarCollection.RcvChannels.ChanParameters.TxRcvPolarization})
    validation_report = add_val_inc(validation_report, 'Error', ...
        'Claimed polarization of image not consistent with what was claimed to have been collected.', ...
        ['SICD_meta.ImageFormation.TxRcvPolarizationProc = ' ...
        SICD_meta.ImageFormation.TxRcvPolarizationProc ...
        ', SICD_meta.RadarCollection.RcvChannels.ChanParameters.TxRcvPolarization: ' ...
        strjoin({SICD_meta.RadarCollection.RcvChannels.ChanParameters.TxRcvPolarization},',')]);    
end

% 2.13 Timeline.IPP
if isfield(SICD_meta,'Timeline') && isfield(SICD_meta.Timeline,'IPP')
    for i = 1:numel(SICD_meta.Timeline.IPP.Set)
        % 2.13.1 IPPStart is IPPPoly evaluated at TStart
        derived_ipp_start = polyval(SICD_meta.Timeline.IPP.Set(i).IPPPoly(end:-1:1), ...
                SICD_meta.Timeline.IPP.Set(i).TStart);
        if abs(derived_ipp_start - SICD_meta.Timeline.IPP.Set(i).IPPStart) > 1
            validation_report = add_val_inc(validation_report, 'Error', ...
                'Timeline.IPP fields not consistent', ...
                ['Timeline.IPP.Set(' num2str(i) ').IPPPoly(Timeline.IPP.Set.TStart): ' num2str(derived_ipp_start), ...
                ', Timeline.IPP.IPPStart: ' num2str(SICD_meta.Timeline.IPP.Set(i).IPPStart)]);
        end
        % 2.13.2 IPPEnd is IPPPoly evaluated at TEnd
        derived_ipp_end = polyval(SICD_meta.Timeline.IPP.Set(i).IPPPoly(end:-1:1), ...
                SICD_meta.Timeline.IPP.Set(i).TEnd);
        if abs(derived_ipp_end - SICD_meta.Timeline.IPP.Set(i).IPPEnd) > 1
            validation_report = add_val_inc(validation_report, 'Error', ...
                'Timeline.IPP fields not consistent', ...
                ['Timeline.IPP.Set(' num2str(i) ').IPPPoly(Timeline.IPP.Set.TEnd): ' num2str(derived_ipp_end), ...
                ', Timeline.IPP.IPPEnd: ' num2str(SICD_meta.Timeline.IPP.Set(i).IPPEnd)]);
        end
    end
    % 2.13.3 All IPP sets should be ordered and consecutive (no skipped IPPs)
    ordered = [SICD_meta.Timeline.IPP.Set.IPPStart] < ...
        [SICD_meta.Timeline.IPP.Set.IPPEnd];
    if numel(SICD_meta.Timeline.IPP.Set)>1
        consecutive = [SICD_meta.Timeline.IPP.Set(1:(end-1)).IPPEnd]+1 == ...
            [SICD_meta.Timeline.IPP.Set(2:end).IPPStart];
    else
        consecutive = true;
    end
    if any(~ordered) || any(~consecutive)
        validation_report = add_val_inc(validation_report, 'Error', ...
            'Timeline.IPP fields not consistent', ...
            'Timeline.IPP.Sets should describe ordered and consecutive IPPs.');
    end
    % 2.13.4 IPP description should match CollectDuration
    if abs((max([SICD_meta.Timeline.IPP.Set.TEnd]) - ...
            min([SICD_meta.Timeline.IPP.Set.TStart])) - ...
            SICD_meta.Timeline.CollectDuration) > .1
        validation_report = add_val_inc(validation_report, 'Error', ...
            'Timeline fields not consistent', ...
            ['Timeline.IPP.Sets span: ' num2str(max([SICD_meta.Timeline.IPP.Set.TEnd]) - ...
            min([SICD_meta.Timeline.IPP.Set.TStart])), ...
            ', Timeline.CollectDuration: ' num2str(SICD_meta.Timeline.CollectDuration)]);
    end
end
% 2.13.5
if (SICD_meta.SCPCOA.SCPTime < 0.25 * SICD_meta.Timeline.CollectDuration) || ...
        (SICD_meta.SCPCOA.SCPTime > 0.75 * SICD_meta.Timeline.CollectDuration)
    validation_report = add_val_inc(validation_report, 'Warning', ...
        'COA time should be roughly in middle of collect.', ...
        ['SICD_meta.SCPCOA.SCPTime: ' num2str(SICD_meta.SCPCOA.SCPTime), ...
        ', Timeline.CollectDuration: ' num2str(SICD_meta.Timeline.CollectDuration)]);
end
% 2.14 Radiometric
if isfield(SICD_meta,'Radiometric')
    % 2.14.1 Relative noise polynomials must have zero as scalar term
    if isfield(SICD_meta.Radiometric,'NoiseLevel') && ...
            strcmpi(SICD_meta.Radiometric.NoiseLevel.NoiseLevelType,'RELATIVE') && ...
            SICD_meta.Radiometric.NoiseLevel.NoisePoly(1)~=0
        validation_report = add_val_inc(validation_report, 'Error', ...
            'Radiometric.NoiseLevel fields not consistent', ...
            ['NoisePoly of type RELATIVE must have a zero scalar term, ' ...
            'since the description is in dB with respective to the SCP.']);
    end
    % Calculate slant plane area
    if isfield(SICD_meta.Grid.Row, 'WgtFunct')
        rng_wght_f=mean(SICD_meta.Grid.Row.WgtFunct.^2) / ...
            (mean(SICD_meta.Grid.Row.WgtFunct).^2);
    else  % If no weight in metadata SICD assumes 1.0
        rng_wght_f=1.0;
    end
    if isfield(SICD_meta.Grid.Col, 'WgtFunct')
        az_wght_f=mean(SICD_meta.Grid.Col.WgtFunct.^2) / ...
            (mean(SICD_meta.Grid.Col.WgtFunct).^2);
    else  % If no weight in metadata SICD assumes 1.0
        az_wght_f=1.0;
    end
    area_sp = (rng_wght_f * az_wght_f) / ...
        (SICD_meta.Grid.Row.ImpRespBW * SICD_meta.Grid.Col.ImpRespBW);
    RAD_POLY_TOL = 1e-5;
    % The RCSSFPoly/SigmaZeroSFPoly/BetaZeroSFPoly/GammaZeroSFPoly fields
    % are all equivalent.
    % 2.14.2
    if isfield(SICD_meta.Radiometric,'RCSSFPoly') && ...
            isfield(SICD_meta.Radiometric,'SigmaZeroSFPoly') && ...
            norm(SICD_meta.Radiometric.RCSSFPoly - ...
            (SICD_meta.Radiometric.SigmaZeroSFPoly*area_sp/...
            cosd(SICD_meta.SCPCOA.SlopeAng))) > RAD_POLY_TOL
        validation_report = add_val_inc(validation_report, 'Error', ...
            'Radiometric.RCSSFPoly/SigmaZeroSFPoly fields not consistent', ...
            ['Radiometric.RCSSFPoly: ' mat2str(SICD_meta.Radiometric.RCSSFPoly) ...
            ', Radiometric.SigmaZeroSFPoly: ' ...
            mat2str(SICD_meta.Radiometric.SigmaZeroSFPoly)]);
    end
    % 2.14.3
    if isfield(SICD_meta.Radiometric,'RCSSFPoly') && ...
            isfield(SICD_meta.Radiometric,'BetaZeroSFPoly') && ...
            norm(SICD_meta.Radiometric.RCSSFPoly - ...
            (SICD_meta.Radiometric.BetaZeroSFPoly*area_sp)) > RAD_POLY_TOL
        validation_report = add_val_inc(validation_report, 'Error', ...
            'Radiometric.RCSSFPoly/BetaZeroSFPoly fields not consistent', ...
            ['Radiometric.RCSSFPoly: ' mat2str(SICD_meta.Radiometric.RCSSFPoly) ...
            ', Radiometric.BetaZeroSFPoly: ' ...
            mat2str(SICD_meta.Radiometric.BetaZeroSFPoly)]);
    end
    % 2.14.4
    if isfield(SICD_meta.Radiometric,'RCSSFPoly') && ...
            isfield(SICD_meta.Radiometric,'GammaZeroSFPoly') && ...
            norm(SICD_meta.Radiometric.RCSSFPoly - ...
            (SICD_meta.Radiometric.GammaZeroSFPoly*area_sp*...
            sind(SICD_meta.SCPCOA.GrazeAng)/...
            cosd(SICD_meta.SCPCOA.SlopeAng))) > RAD_POLY_TOL
        validation_report = add_val_inc(validation_report, 'Error', ...
            'Radiometric.RCSSFPoly/GammaZeroSFPoly fields not consistent', ...
            ['Radiometric.RCSSFPoly: ' mat2str(SICD_meta.Radiometric.RCSSFPoly) ...
            ', Radiometric.GammaZeroSFPoly: ' ...
            mat2str(SICD_meta.Radiometric.GammaZeroSFPoly)]);
    end
    % 2.14.5
    if isfield(SICD_meta.Radiometric,'SigmaZeroSFPoly') && ...
            isfield(SICD_meta.Radiometric,'BetaZeroSFPoly') && ...
            norm(SICD_meta.Radiometric.SigmaZeroSFPoly - ...
            (SICD_meta.Radiometric.BetaZeroSFPoly*...
            cosd(SICD_meta.SCPCOA.SlopeAng))) > RAD_POLY_TOL
        validation_report = add_val_inc(validation_report, 'Error', ...
            'Radiometric.SigmaZeroSFPoly/BetaZeroSFPoly fields not consistent', ...
            ['Radiometric.SigmaZeroSFPoly: ' mat2str(SICD_meta.Radiometric.SigmaZeroSFPoly) ...
            ', Radiometric.BetaZeroSFPoly: ' ...
            mat2str(SICD_meta.Radiometric.BetaZeroSFPoly)]);
    end
    % 2.14.6
    if isfield(SICD_meta.Radiometric,'SigmaZeroSFPoly') && ...
            isfield(SICD_meta.Radiometric,'GammaZeroSFPoly') && ...
            norm(SICD_meta.Radiometric.SigmaZeroSFPoly - ...
            (SICD_meta.Radiometric.GammaZeroSFPoly*...
            sind(SICD_meta.SCPCOA.GrazeAng))) > RAD_POLY_TOL
        validation_report = add_val_inc(validation_report, 'Error', ...
            'Radiometric.SigmaZeroSFPoly/GammaZeroSFPoly fields not consistent', ...
            ['Radiometric.SigmaZeroSFPoly: ' mat2str(SICD_meta.Radiometric.SigmaZeroSFPoly) ...
            ', Radiometric.GammaZeroSFPoly: ' ...
            mat2str(SICD_meta.Radiometric.GammaZeroSFPoly)]);
    end
    % 2.14.7
    if isfield(SICD_meta.Radiometric,'BetaZeroSFPoly') && ...
            isfield(SICD_meta.Radiometric,'GammaZeroSFPoly') && ...
            norm(SICD_meta.Radiometric.BetaZeroSFPoly - ...
            (SICD_meta.Radiometric.GammaZeroSFPoly*...
            sind(SICD_meta.SCPCOA.GrazeAng)/...
            cosd(SICD_meta.SCPCOA.SlopeAng))) > RAD_POLY_TOL
        validation_report = add_val_inc(validation_report, 'Error', ...
            'Radiometric.BetaZeroSFPoly/GammaZeroSFPoly fields not consistent', ...
            ['Radiometric.BetaZeroSFPoly: ' mat2str(SICD_meta.Radiometric.BetaZeroSFPoly) ...
            ', Radiometric.GammaZeroSFPoly: ' ...
            mat2str(SICD_meta.Radiometric.GammaZeroSFPoly)]);
    end
end

%% 3 Check for consistency between SICD metadata and pixel data through interactive visual checks
if do_interactive
    % 3.1 View spatial variance of frequency support and overlay
    % predicted frequency support from
    % Row/Col.DeltaKCOAPoly/ImpRespBW/DeltaK1/DeltaK2
    % This is a quick check to manually check spatially variant frequency
    % support.  An option is provided later for a more complete full-image
    % check (slower since it has to process the entire dataset).
    if ~fs_vis_test(sicd_filename)
        validation_report = add_val_inc(validation_report, 'Error', ...
            'SICD metadata incorrectly describe spatial frequency bounds.', ...
            'DeltaKCOAPoly values manually reviewed.');
    end
    % 3.2 Use SICD sensor model to ground project a low resolution version
    % of the amplitude image.  Create KML/KMZ to overlay on GoogleEarth to
    % assure orientation and rough geolocation of the image is correct.
    k = kml(SICD_meta.CollectionInfo.CoreName); %create an kmltoolbox object
    k.filename = tempname();
    k.zip = true; % We want to wrap the overlay image with KML into a single KMZ
    add_sar_2kml(k,sicd_filename,'overlay_decimate','max',...
        'overlay_max_size',1000,'overlay_filename','overlay.png');
    k.save();
    delete(k.includeFiles{:}); % Overlay file
    try
        winopen(k.filename);
        display(['KML file temporarily saved at ' k.filename '.']);
    catch
        display(['Unable to automatically open Google Earth.  Please ' ...
            'manually open file temporarily saved at ' k.filename ...
            ' in Google Earth to check imagery orientation and rough geolocation.']);
    end
    button = questdlg('Does data overlayed in Google Earth match background imagery?','Amplitude Image Verification','Yes','No','Yes');
    if strcmp(button, 'No')
        validation_report = add_val_inc(validation_report, 'Error', ...
            'Image orientation or geolocation incorrect.', ...
            'Amplitude image manually reviewed.');
    end
    delete(k.filename); % Remove KMZ file when done with it
    % 3.3 Visualize travelling glint to check FFT sign
    button = questdlg('Does this dataset have a travelling glint or other feature that would allow FFT sign verication?','FFT Sign Verification','Yes','No','Yes');
    if strcmp(button, 'Yes')
        % The ApertureTool is the best for this.  Requires a collect with a
        % travelling glint (or perhaps a shadow for high-resolution
        % collections) for verification though.  ApertureTool will process
        % with respect to the FFT sign given in the file.  If a slow-time
        % subaperture animation shows travelling glint or shadow moving in
        % the expected direction, with respect to SideOfTrack, then the FFT
        % sign is populated correctly.
        % Another way to check this would be in fast-time, instead of
        % slow-time.  This would require a visible (not-inscribed) polar
        % annulus or a known narrowband emitter in the data.
        aptool_fig = ApertureTool('filename',sicd_filename);
        uiwait(aptool_fig);
        button = questdlg('Did subaperture signatures behave correctly for flight direction?','FFT Sign Verification','Yes','No','Yes');
        if strcmp(button, 'No')
            validation_report = add_val_inc(validation_report, 'Error', ...
                'FFT sign incorrect.', ...
                ['Grid.Row/Col.Sgn: ' num2str(SICD_meta.Grid.Col.Sgn)]); % Row and Col should be the same
        end
    end
    button = questdlg('Do you want to do full-image frequency support analysis?  This can take a while on very large datasets.',...
        'Spatial Frequency Support Verification','Yes','No','Yes');
    if strcmp(button, 'Yes')
        % 3.4 Visualize image domain spatial frequency description
        FS_VIS_SIZE = 1000; % We want to view spatial in an image about this size
        % 3.4.1 DeltaK1/DeltaK2 check
        fft_filename = tempname();
        fftfile(sicd_filename, fft_filename, 'azlimits', 'full', 'rnglimits', 'full'); % Compute full FFT of entire file
        decimation_level = max(1, ceil(double([SICD_meta.ImageData.NumCols SICD_meta.ImageData.NumRows])./FS_VIS_SIZE));
        full_fft = read_complex_data(fft_filename, [], [], decimation_level, 'max');
        delete(fft_filename); % Don't need full file anymore
        % Note that fftfile.m does its computations with the FFT sign
        % intentionally inverted. This was done so that a typical MATLAB
        % display (which typically has YDir='reverse' for image axes) would
        % result in RF frequency increasing upward and polar angles moving in
        % the clockwise direction as you move toward the right.  However the
        % convention in SICD is RF frequency increases with increasing Krow
        % index (down in a traditional MATLAB display) and polar angle moves in
        % the counterclockwise direction with increasing Kcol index (which is
        % consistent with the polar angle convetion).  So we flip the result
        % here so that the SICD conventions follow the MATLAB indexing.
        full_fft = rot90(full_fft,2);
        full_fig = figure('Name','Full FFT');
        full_im = imshow(logremap(full_fft.'));
        full_ax = get(full_im,'Parent');
        set(full_ax,'YDir','normal'); % So that Krow increases up
        xlabel(full_ax, 'Kcol');
        ylabel(full_ax, 'Krow');
        title(full_ax, 'Full FFT');
        % The following are computed in the units of the MATLAB axes, which are
        % based off of pixels in the displayed image.
        krow_bounds = [round((size(full_fft,2)-1) * ...
            (0.5 + SICD_meta.Grid.Row.SS*SICD_meta.Grid.Row.DeltaK1)) + 1, ...
            round((size(full_fft,2)-1) * ...
            (0.5 + SICD_meta.Grid.Row.SS*SICD_meta.Grid.Row.DeltaK2)) + 1];
        kcol_bounds = [round((size(full_fft,1)-1) * ...
            (0.5 + SICD_meta.Grid.Col.SS*SICD_meta.Grid.Col.DeltaK1)) + 1, ...
            round((size(full_fft,1)-1) * ...
            (0.5 + SICD_meta.Grid.Col.SS*SICD_meta.Grid.Col.DeltaK2)) + 1];
        hold(full_ax,'on');
        plot(full_ax, xlim(full_ax), krow_bounds(1)*[1 1], '--r');
        plot(full_ax, xlim(full_ax), krow_bounds(2)*[1 1], '--r');
        plot(full_ax, kcol_bounds(1)*[1 1], ylim(full_ax), '--r');
        plot(full_ax, kcol_bounds(2)*[1 1], ylim(full_ax), '--r');
        set(full_ax, 'Visible', 'on');
        try % Don't crash for DeltaK1=DeltaK2 case
            set(full_ax, 'XTick', [min(kcol_bounds) mean(xlim(full_ax)) max(kcol_bounds)]);
            set(full_ax, 'XTickLabel', {'\DeltaKcol1', '0', '\DeltaKcol2'});
            set(full_ax, 'YTick', [min(krow_bounds) mean(ylim(full_ax)) max(krow_bounds)]);
            set(full_ax, 'YTickLabel', {'\DeltaKrow1', '0', '\DeltaKrow2'});
        catch
            set(full_ax, 'Visible', 'off');
        end
        button = questdlg('Does displayed data match Krow and KCol bounds shown in red?','Spatial Frequency Support Verification','Yes','No','Yes');
        if strcmp(button, 'No')
            validation_report = add_val_inc(validation_report, 'Error', ...
                'SICD metadata incorrectly describe spatial frequency bounds.', ...
                'DeltaK1/DeltaK2 values manually reviewed.');
        end
        % 3.4.2 Col.DeltaKCOAPoly/ImpRespBW
        col_bw_min = round((size(full_fft,1)-1) * ...
            (0.5 - SICD_meta.Grid.Col.SS*SICD_meta.Grid.Col.ImpRespBW/2)) + 1;
        col_bw_max = round((size(full_fft,1)-1) * ...
            (0.5 + SICD_meta.Grid.Col.SS*SICD_meta.Grid.Col.ImpRespBW/2)) + 1;
        if isfield(SICD_meta.Grid.Col, 'DeltaKCOAPoly') && any(SICD_meta.Grid.Col.DeltaKCOAPoly(:))
            col_filename = tempname();
            deskewfile(sicd_filename, 'output_filename', col_filename, 'dim', 1);
            col_fft_filename = tempname();
            fftfile(col_filename, col_fft_filename, 'azlimits', 'full', ...
                'rnglimits', 'full', 'sgn', SICD_meta.Grid.Col.Sgn);
            col_fft = read_complex_data(col_fft_filename, [], [], decimation_level, 'max');
            col_fft = rot90(col_fft,2); % See comments in test 3.1.1 for why
            figure('Name','KCol aligned FFT');
            col_im = imshow(logremap(col_fft.'));
            col_ax = get(col_im,'Parent');
            set(col_ax,'YDir','normal'); % So that Krow increases up
            xlabel(col_ax, 'Kcol');
            ylabel(col_ax, 'Krow');
            title(col_ax, 'FFT of KCol aligned data');
            hold(col_ax,'on');
            plot(col_ax, [col_bw_min col_bw_min], ylim(col_ax), '--b');
            plot(col_ax, [col_bw_max col_bw_max], ylim(col_ax), '--b');
            set(col_ax, 'Visible', 'on');
            set(col_ax, 'YTick', []);
            try % Don't crash for ImpRespBW=0 case
                set(col_ax, 'XTick', [col_bw_min mean(xlim(col_ax)) col_bw_max]);
                set(col_ax, 'XTickLabel', {'-Kcol_{IRBW}/2', '0', 'Kcol_{IRBW}/2'});
            end
            delete(col_filename);
            delete(col_fft_filename);
            col_wgt = sum(abs(col_fft),2);
        else
            figure(full_fig);
            plot(full_ax, [col_bw_min col_bw_min], ylim(full_ax), '--b');
            plot(full_ax, [col_bw_max col_bw_max], ylim(full_ax), '--b');
            set(full_ax, 'YTick', []);
            try % Don't crash for ImpRespBW=0 case
                set(full_ax, 'XTick', [col_bw_min mean(xlim(full_ax)) col_bw_max]);
                set(full_ax, 'XTickLabel', {'-Kcol_{IRBW}/2', '0', 'Kcol_{IRBW}/2'});
            end
            col_wgt = sum(abs(full_fft),2);
        end
        button = questdlg('Does displayed data match Col.ImpRespBW bounds shown in blue?','Spatial Frequency Support Verification','Yes','No','Yes');
        if strcmp(button, 'No')
            validation_report = add_val_inc(validation_report, 'Error', ...
                'SICD metadata incorrectly describe spatial frequency bounds.', ...
                'Col.DeltaKCOAPoly/ImpRespBW values manually reviewed.');
        end
        % 3.4.3 Row.DeltaKCOAPoly/ImpRespBW
        row_bw_min = round((size(full_fft,2)-1) * ...
            (0.5 - SICD_meta.Grid.Row.SS*SICD_meta.Grid.Row.ImpRespBW/2)) + 1;
        row_bw_max = round((size(full_fft,2)-1) * ...
            (0.5 + SICD_meta.Grid.Row.SS*SICD_meta.Grid.Row.ImpRespBW/2)) + 1;
        if isfield(SICD_meta.Grid.Row, 'DeltaKCOAPoly') && any(SICD_meta.Grid.Row.DeltaKCOAPoly(:))
            row_filename = tempname();
            deskewfile(sicd_filename, 'output_filename', row_filename, 'dim', 2);
            row_fft_filename = tempname();
            fftfile(row_filename, row_fft_filename, 'azlimits', 'full', ...
                'rnglimits', 'full', 'sgn', SICD_meta.Grid.Row.Sgn);
            row_fft = read_complex_data(row_fft_filename, [], [], decimation_level, 'max');
            row_fft = rot90(row_fft,2); % See comments in test 3.1.1 for why
            figure('Name','KRow aligned FFT');
            row_im = imshow(logremap(row_fft.'));
            row_ax = get(row_im,'Parent');
            set(row_ax,'YDir','normal'); % So that Krow increases up
            xlabel(row_ax, 'Kcol');
            ylabel(row_ax, 'Krow');
            title(row_ax, 'FFT of KRow aligned data');
            hold(row_ax,'on');
            plot(row_ax, xlim(full_ax), [row_bw_min row_bw_min], '--g');
            plot(row_ax, xlim(full_ax), [row_bw_max row_bw_max], '--g');
            set(row_ax, 'Visible', 'on');
            set(row_ax, 'XTick', []);
            try % Don't crash for ImpRespBW=0 case
                set(row_ax, 'YTick', [row_bw_min mean(ylim(row_ax)) row_bw_max]);
                set(row_ax, 'YTickLabel', {'-Krow_{IRBW}/2', '0', 'Krow_{IRBW}/2'});
            end
            delete(row_filename);
            delete(row_fft_filename);
            row_wgt = sum(abs(row_fft));
        else
            figure(full_fig);
            plot(full_ax, xlim(full_ax), [row_bw_min row_bw_min], '--g');
            plot(full_ax, xlim(full_ax), [row_bw_max row_bw_max], '--g');
            set(full_ax, 'XTick', []);
            try % Don't crash for ImpRespBW=0 case
                set(full_ax, 'YTick', [row_bw_min mean(ylim(full_ax)) row_bw_max]);
                set(full_ax, 'YTickLabel', {'-Krow_{IRBW}/2', '0', 'Krow_{IRBW}/2'});
            end
            row_wgt = sum(abs(full_fft));
        end
        button = questdlg('Does displayed data match Row.ImpRespBW bounds shown in green?','Spatial Frequency Support Verification','Yes','No','Yes');
        if strcmp(button, 'No')
            validation_report = add_val_inc(validation_report, 'Error', ...
                'SICD metadata incorrectly describe spatial frequency bounds.', ...
                'Row.DeltaKCOAPoly/ImpRespBW values manually reviewed.');
        end
        % 3.5 Measure and display weighting across image and compared to
        % claimed weighting.
        % 3.5.1 Row
        figure('Name', 'Row Weighting');
        row_wgt_ax = gca();
        scale_factor = computeWgtSF(row_wgt, ...
            SICD_meta.Grid.Row.ImpRespBW*SICD_meta.Grid.Row.SS);
        plot(row_wgt_ax, ...
            linspace(-0.5/SICD_meta.Grid.Row.SS, 0.5/SICD_meta.Grid.Row.SS, numel(row_wgt)), ...
            row_wgt*scale_factor, '-b');
        xlabel(row_wgt_ax, 'Krow (cycles/meter)');
        hold(row_wgt_ax,'on');
        legend_str = {'From complex data'};
        if exist('rowFun','var')
            if ~isempty(rowFun)
                rowFunSamples = rowFun(DEFAULT_WGT_SIZE);
            else % Uniform weighting
                rowFunSamples = ones(2,1);
            end
            plot(row_wgt_ax, ...
                linspace(-SICD_meta.Grid.Row.ImpRespBW/2, ...
                SICD_meta.Grid.Row.ImpRespBW/2, ...
                numel(rowFunSamples)), rowFunSamples, '-r');
            legend_str{end+1} = 'Grid.Row.WgtType';
        end
        if isfield(SICD_meta.Grid.Row,'WgtFunct')
            plot(row_wgt_ax, ...
                linspace(-SICD_meta.Grid.Row.ImpRespBW/2, ...
                SICD_meta.Grid.Row.ImpRespBW/2, ...
                numel(SICD_meta.Grid.Row.WgtFunct)), ...
                SICD_meta.Grid.Row.WgtFunct, '--g');
            legend_str{end+1} = 'Grid.Row.WgtFunct';
        end
        legend(legend_str{:});
        button = questdlg('Do amplitude weightings derived from pixel data and XML metadata agree?','Weighting Verification','Yes','No','Yes');
        if strcmp(button, 'No')
            validation_report = add_val_inc(validation_report, 'Error', ...
                'SICD metadata incorrectly describes amplitude weighting seen in data.', ...
                'Row.WgtType/WgtFunct values manually reviewed.');
        end
        % 3.5.2 Column
        figure('Name', 'Column Weighting');
        col_wgt_ax = gca();
        scale_factor = computeWgtSF(col_wgt, ...
            SICD_meta.Grid.Col.ImpRespBW*SICD_meta.Grid.Col.SS);
        plot(col_wgt_ax, ...
            linspace(-0.5/SICD_meta.Grid.Col.SS, 0.5/SICD_meta.Grid.Col.SS, numel(col_wgt)), ...
            col_wgt*scale_factor, '-b');
        xlabel(col_wgt_ax, 'Kcol (cycles/meter)');
        hold(col_wgt_ax,'on');
        legend_str = {'From complex data'};
        if exist('colFun','var')
            if ~isempty(colFun)
                colFunSamples = colFun(DEFAULT_WGT_SIZE);
            else % Uniform weighting
                colFunSamples = ones(2,1);
            end
            plot(col_wgt_ax, ...
                linspace(-SICD_meta.Grid.Col.ImpRespBW/2, ...
                SICD_meta.Grid.Col.ImpRespBW/2, ...
                numel(colFunSamples)), colFunSamples, '-r');
            legend_str{end+1} = 'Grid.Col.WgtType';
        end
        if isfield(SICD_meta.Grid.Col,'WgtFunct')
            plot(col_wgt_ax, ...
                linspace(-SICD_meta.Grid.Col.ImpRespBW/2, ...
                SICD_meta.Grid.Col.ImpRespBW/2, ...
                numel(SICD_meta.Grid.Col.WgtFunct)), ...
                SICD_meta.Grid.Col.WgtFunct, '--g');
            legend_str{end+1} = 'Grid.Col.WgtFunct';
        end
        legend(legend_str{:});
        button = questdlg('Do amplitude weightings derived from pixel data and XML metadata agree?','Weighting Verification','Yes','No','Yes');
        if strcmp(button, 'No')
            validation_report = add_val_inc(validation_report, 'Error', ...
                'SICD metadata incorrectly describes amplitude weighting seen in data.', ...
                'Col.WgtType/WgtFunct values manually reviewed.');
        end
    end
    % 3.6 Compare measured noise to predicted noise from metadata
    if isfield(SICD_meta, 'Radiometric') && ...
            isfield(SICD_meta.Radiometric,'NoiseLevel') && ...
            strcmpi(SICD_meta.Radiometric.NoiseLevel.NoiseLevelType,'ABSOLUTE')
        button = questdlg('Does this dataset have low return areas that could be selected for noise level verication?','Noise Verification','Yes','No','Yes');
        if strcmp(button, 'Yes')
            test_noise(sicd_filename);
            uiwait(gcf);
            button = questdlg('Was measured noise sufficiently close to predicted noise?','Noise Verification','Yes','No','Yes');
            if strcmp(button, 'No')
                validation_report = add_val_inc(validation_report, 'Error', ...
                    'Noise level description incorrect.', ...
                    'Radiometric.NoiseLevel was manually reviewed.');
            end
        end
    end
end

%% 4 Best practices check
% 4.1
% Extremely long WgtFunct are clumsy and unecessary.  They can result in
% slow XML parsing, and they provide no additional information since
% weighting functions are always quite smooth.
if (isfield(SICD_meta.Grid.Row, 'WgtFunct') && ...
        numel(SICD_meta.Grid.Row.WgtFunct)>1024) || ...
        (isfield(SICD_meta.Grid.Row, 'WgtFunct') && ...
        numel(SICD_meta.Grid.Col.WgtFunct)>1024)
    validation_report = add_val_inc(validation_report, 'Style', ...
        'Weighting functions longer than necessary.', ...
        'Recommended length of Grid.Row/Col.WgtFunct is 512.');
end
% 4.2 Optional fields that are frequently useful:
% 4.2.1 ErrorStatistics
if ~isfield(SICD_meta, 'ErrorStatistics')
    validation_report = add_val_inc(validation_report, 'Optional', ...
        'No ErrorStatistics field provided.', ...
        'No error propogation through a rigorous sensor model will be available.');
end
% 4.2.2 Radiometric
if ~isfield(SICD_meta, 'Radiometric') || ~any([...
        isfield(SICD_meta.Radiometric,'RCSSFPoly') ... % We generally use this one, but we belive them to all be equivalent
        isfield(SICD_meta.Radiometric,'SigmaZeroSFPoly') ...
        isfield(SICD_meta.Radiometric,'BetaZeroSFPoly') ...
        isfield(SICD_meta.Radiometric,'GammaZeroSFPoly')])
    validation_report = add_val_inc(validation_report, 'Optional', ...
        'No radiometric calibration fields provided.', ...
        'No absolute RCS measurements will be available in this data.');
end
% 4.2.3 Noise
if ~isfield(SICD_meta, 'Radiometric') || ...
        ~isfield(SICD_meta.Radiometric, 'NoiseLevel') || ...
        ~strcmp(SICD_meta.Radiometric.NoiseLevel.NoiseLevelType, 'ABSOLUTE')
    validation_report = add_val_inc(validation_report, 'Optional', ...
        'No absolute noise values provided.', ...
        'Absolute noise values useful for a number of SAR analysis techniques including some interferometric calculations.');
end
% 4.2.4 Timeline.IPP
if ~isfield(SICD_meta, 'Timeline') || ...
        ~isfield(SICD_meta.Timeline, 'IPP')
    validation_report = add_val_inc(validation_report, 'Optional', ...
        'No output PRF/PRI info provided.', ...
        'Timeline.IPP would allow for better analysis of ambiguities.');
end
% 4.2.5 RadarCollection.Area.Plane
if ~isfield(SICD_meta, 'RadarCollection') || ...
        ~isfield(SICD_meta.RadarCollection, 'Area') || ...
        ~isfield(SICD_meta.RadarCollection.Area, 'Plane') 
    validation_report = add_val_inc(validation_report, 'Optional', ...
        'No output Area/Plane provided.', ...
        'Some tools prefer using a predetermined output plane layout for consistent output products.');
end
% 4.2.6 ValidData
if ~isfield(SICD_meta.ImageData, 'ValidData')
    validation_report = add_val_inc(validation_report, 'Optional', ...
        'No ValidData area provided.', ...
        ['If there is any chance for pixels in the image that are ' ...
        'zeroed out or do not contain valid image information, then this ' ...
        'field should be filled out.  If it is assured that all pixels ' ...
        'are valid, it is still recommended to populate this explicitly.']);
end
% 4.2.7 RefFreqIndex
if isfield(SICD_meta.RadarCollection,'RefFreqIndex')
    validation_report = add_val_inc(validation_report, 'Optional', ...
        'Reference frequency is being used.', ...
        ['All RF frequency values are expressed as offsets from a ' ...
        'reference frequency.  A number of validation tests could not ' ...
        'be performed because of this.']);
end
% RadarCollection.Waveform is usually populated, and we could check for it
% as well, but it seems silly to throw a warning for a field not used in
% exploitation.  Likewise the Antenna field is equally as useless for
% exploitation.

catch ME
    validation_report = add_val_inc(validation_report, 'Error', ...
        'SICD validation crashed.', getReport(ME,'extended','hyperlinks','off'));
end

end

% Add validation incident
function [ validation_report ] = add_val_inc( validation_report, level, message, details )
    validation_report{end+1,1} = level;
    validation_report{end,2} = message;
    validation_report{end,3} = details;
end

function ImpRespWid = computeImpRespWid(WgtFunct, ImpRespBW)
    OVERSAMPLE = 1024;
    if isempty(WgtFunct) % Way to denote uniform weighting
        broadening_factor = 2 * fzero(@(x) (sin(pi*x)/(pi*x)) - (1/sqrt(2)), .1); % 0.886
    else
        imp_resp = abs(fft(WgtFunct, round(numel(WgtFunct)*OVERSAMPLE))); % Oversampled response function
        imp_resp = imp_resp/sum(WgtFunct); % Normalize to unit peak
        ind = find(imp_resp<1/sqrt(2),1,'first')+[-1 -0]; % Samples surrounding half-power point
        ind = interp1(imp_resp(ind), ind, 1/sqrt(2)); % Linear interpolation to solve for half-power point
        broadening_factor = 2*(ind - 1)/OVERSAMPLE;
    end
    ImpRespWid = broadening_factor/ImpRespBW;
end

% We try to compute a scale factor that will appropriately scale a
% data-derived weighting (of arbitrary amplitude) to a peak of 1.  We do
% this slightly differently for low-pass weightings vs. uniform
% weightings.
function scale_factor = computeWgtSF(profile1d, freq_support_percent)
    profile1d = profile1d(:); % Assure consistent orientation
    freq_support_percent = max(min(freq_support_percent,1),0); % Assure reasonable value
    % We create a very low pass version of data-derived weighting.
    % Typically we don't need that many components to approximate a
    % weighting.  Many weightings (like Hamming/Hanning) are described with
    % only 2 (non-negative) FFT elements.  The number we use here is
    % totally arbitrary with no formal justification.
    LP_FFT_ELEMS = 7; % Number of non-negative FFT sinusoidal components to use for low-pass estimate
    trim = ceil(numel(profile1d) * (1-freq_support_percent) / 2);
    trimmed_data = profile1d((trim+1):(end-trim-1));
    sinusoidal_fit = fft(trimmed_data);
    sinusoidal_fit((LP_FFT_ELEMS+1):(end-LP_FFT_ELEMS+1)) = 0;
    sinusoidal_fit = ifft(sinusoidal_fit);
    sinusoidal_peak = max(sinusoidal_fit);
    % We assume any non-uniform weighting will peak in the middle and taper
    % down by some signifcant amount by the edges.
    if (sinusoidal_peak*.8)> max(sinusoidal_fit([1 end])) % Likely a 
        scale_factor = 1/max(sinusoidal_fit); % Make max of smoothed data 1
    else % Likely a uniform weighting
        scale_factor = 1/mean(trimmed_data); % Make mean of data 1
    end
end

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////